oxlith_manuscript

theme_set(theme_classic()+
          theme(text = element_text(family = "Helvetica", size = 10)))
orangeNavyPalette = c("orange","steelblue4")

df_hlm_movement <- read.csv("../data/movement_data_HLM_-14toEnd_allparticipants.csv")
names(df_hlm_movement)[6:13] <- paste0("movement_",names(df_hlm_movement)[6:13])
df_hlm_time <- read.csv("../data/time_data_HLM_-14toEnd_allparticipants.csv")
names(df_hlm_time)[6:13] <- paste0("time_",names(df_hlm_time)[6:13])
df_hlm_mood <- read.csv("../data/mood_data_HLM_-14toEnd_allparticipants.csv")

df_hlm <- Reduce(merge, list(df_hlm_movement,df_hlm_time,df_hlm_mood))

df_demo <- read.csv("../data/demo_info.csv")%>%
  mutate(Allocation = if_else(Allocation == "B","Lithium", "Control"))
  # filter(ID %in% unique(df_hlm$ID))

participant_num <- length(unique(df_hlm$ID))

df_hlm$Allocation[df_hlm$Allocation == "A"] <- "control"
df_hlm$Allocation[df_hlm$Allocation == "B"] <- "lithium"

# mood parameters from 0~4 to 1~5
df_hlm$pos_mean <- df_hlm$pos_mean+5
df_hlm$neg_mean <- df_hlm$neg_mean+5

# mood parameters need log transformation
df_hlm$posLog <- log(df_hlm$pos_mean+1)
df_hlm$negLog <- log(df_hlm$neg_mean+1)

#compute relative amplitude
df_hlm <- df_hlm%>%
  rowwise()%>%
  mutate(rel_amplitude  = (movement_M10_mean - movement_L5_mean)/(movement_M10_mean + movement_L5_mean))

#add season into df
df_hlm <- merge(df_hlm,df_demo[,c("ID","Season","Diagnosis")], by = "ID")

demo descriptive stats

#reorganize ethnicity
df_demo <- df_demo%>%
  mutate(Ethnicity = if_else(Ethnicity == "", "not reported", if_else(Ethnicity == "other brazilian" | Ethnicity == "Other croatian british","other",Ethnicity)))

df_control <- df_demo%>%filter(Allocation == "Control")
df_lithium <- df_demo%>%filter(Allocation == "Lithium")

table(df_demo[,c("Gender","Allocation")])

##       Allocation
## Gender Control Lithium
##      F       9      11
##      M       7       8

chisq.test(table(df_demo[,c("Gender","Allocation")]))

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(df_demo[, c("Gender", "Allocation")])
## X-squared = 1.1504e-31, df = 1, p-value = 1

#ethnicity
table(df_demo[,c("Ethnicity","Allocation")])

##                 Allocation
## Ethnicity        Control Lithium
##   Asian                1       0
##   Black american       0       1
##   Hispanic             0       1
##   mixed                1       1
##   not reported         3       1
##   other                1       1
##   wb                   9      13
##   wo                   1       1

chisq.test(table(df_demo[,c("Ethnicity","Allocation")]))

## Warning in chisq.test(table(df_demo[, c("Ethnicity", "Allocation")])):
## Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  table(df_demo[, c("Ethnicity", "Allocation")])
## X-squared = 4.5032, df = 7, p-value = 0.7203

#age
df_demo%>%
  group_by(Allocation)%>%
  summarise(mean_Age = round(mean(Age),2),
            sd_Age = round(sd(Age),2))

## # A tibble: 2 × 3
##   Allocation mean_Age sd_Age
##   <chr>         <dbl>  <dbl>
## 1 Control        34.7  13.9 
## 2 Lithium        28.4   9.86

t.test(df_control$Age,df_lithium$Age)

## 
##  Welch Two Sample t-test
## 
## data:  df_control$Age and df_lithium$Age
## t = 1.5269, df = 26.485, p-value = 0.1386
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.180022 14.818180
## sample estimates:
## mean of x mean of y 
##  34.68750  28.36842

#BMI
df_demo%>%
  group_by(Allocation)%>%
  summarise(mean_BMI = round(mean(BMI),2),
            sd_BMI = round(sd(BMI),2))

## # A tibble: 2 × 3
##   Allocation mean_BMI sd_BMI
##   <chr>         <dbl>  <dbl>
## 1 Control        26.5   4.95
## 2 Lithium        25.4   6.28

t.test(df_control$BMI,df_lithium$BMI)

## 
##  Welch Two Sample t-test
## 
## data:  df_control$BMI and df_lithium$BMI
## t = 0.57297, df = 32.884, p-value = 0.5706
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.777079  4.954053
## sample estimates:
## mean of x mean of y 
##  26.49375  25.40526

#mood
df_mood <- df_hlm%>%
  filter(day >= -14 & day <= 28)%>%
  group_by(ID,Allocation)%>%
  summarise(pos_mean = round(mean(pos_mean,na.rm = TRUE),2),
            neg_mean = round(mean(neg_mean,na.rm = TRUE),2))

## `summarise()` has grouped output by 'ID'. You can override using the `.groups`
## argument.

df_mood%>%
  group_by(Allocation)%>%
  summarise(mean_pos = round(mean(pos_mean),2),
            sd_po = round(sd(pos_mean),2),
            mean_neg = round(mean(neg_mean),2),
            sd_neg = round(sd(neg_mean),2))

## # A tibble: 2 × 5
##   Allocation mean_pos sd_po mean_neg sd_neg
##   <chr>         <dbl> <dbl>    <dbl>  <dbl>
## 1 control        10.4  3.01    10.7    4.18
## 2 lithium        10.7  3.5      8.71   3.25

t.test(df_mood$pos_mean[df_mood$Allocation == "lithium"],df_mood$pos_mean[df_mood$Allocation == "control"])

## 
##  Welch Two Sample t-test
## 
## data:  df_mood$pos_mean[df_mood$Allocation == "lithium"] and df_mood$pos_mean[df_mood$Allocation == "control"]
## t = 0.24204, df = 31.708, p-value = 0.8103
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.006914  2.547967
## sample estimates:
## mean of x mean of y 
##  10.71053  10.44000

t.test(df_mood$neg_mean[df_mood$Allocation == "lithium"],df_mood$neg_mean[df_mood$Allocation == "control"])

## 
##  Welch Two Sample t-test
## 
## data:  df_mood$neg_mean[df_mood$Allocation == "lithium"] and df_mood$neg_mean[df_mood$Allocation == "control"]
## t = -1.4999, df = 25.943, p-value = 0.1457
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -4.6624583  0.7288442
## sample estimates:
## mean of x mean of y 
##  8.710526 10.677333

#recording days
df_demo%>%
  group_by(Allocation)%>%
  summarise(median_total.recording.days = median(total.recording.days, na.rm = TRUE))

## # A tibble: 2 × 2
##   Allocation median_total.recording.days
##   <chr>                            <dbl>
## 1 Control                             50
## 2 Lithium                             50

table(df_demo[,c("total.recording.days","Allocation")])

##                     Allocation
## total.recording.days Control Lithium
##                   46       2       0
##                   47       2       1
##                   48       0       4
##                   49       2       2
##                   50       7      10
##                   52       0       1

t.test(df_lithium$total.recording.days,df_control$total.recording.days)

## 
##  Welch Two Sample t-test
## 
## data:  df_lithium$total.recording.days and df_control$total.recording.days
## t = 1.1578, df = 20.801, p-value = 0.2601
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.4940028  1.7333190
## sample estimates:
## mean of x mean of y 
##  49.38889  48.76923

#missing days
df_demo%>%
  group_by(Allocation)%>%
  summarise(median_Missing.days = median(Missing.days, na.rm = TRUE))

## # A tibble: 2 × 2
##   Allocation median_Missing.days
##   <chr>                    <dbl>
## 1 Control                   17  
## 2 Lithium                   13.5

table(df_demo[,c("Missing.days","Allocation")])

##             Allocation
## Missing.days Control Lithium
##           0        1       3
##           3        1       0
##           4        1       2
##           7        0       1
##           9        0       1
##           11       1       1
##           12       1       1
##           14       1       0
##           15       0       1
##           17       2       0
##           18       0       1
##           19       1       1
##           22       2       0
##           24       1       0
##           25       0       1
##           27       0       1
##           31       0       1
##           32       0       1
##           34       1       0
##           35       0       2

t.test(df_lithium$Missing.days,df_control$Missing.days)

## 
##  Welch Two Sample t-test
## 
## data:  df_lithium$Missing.days and df_control$Missing.days
## t = 0.11917, df = 28.888, p-value = 0.906
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -7.598748  8.538919
## sample estimates:
## mean of x mean of y 
##  15.77778  15.30769

#BP subtype
table(df_demo[,c("Diagnosis","Allocation")])

##          Allocation
## Diagnosis Control Lithium
##       1         4       3
##       2        11      16
##       NOS       1       0

chisq.test(table(df_demo[,c("Diagnosis","Allocation")]))

## Warning in chisq.test(table(df_demo[, c("Diagnosis", "Allocation")])):
## Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  table(df_demo[, c("Diagnosis", "Allocation")])
## X-squared = 1.825, df = 2, p-value = 0.4015

#season
table(df_demo[,c("Season","Allocation")])

##         Allocation
## Season   Control Lithium
##   Autumn       4       4
##   Spring       3       3
##   Summer       5       3
##   Winter       4       9

chisq.test(table(df_demo[,c("Season","Allocation")]))

## Warning in chisq.test(table(df_demo[, c("Season", "Allocation")])): Chi-squared
## approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  table(df_demo[, c("Season", "Allocation")])
## X-squared = 2.182, df = 3, p-value = 0.5355

data inspection

data availability

participant_num

## [1] 34

data_path <- "../data/"
df_hlm_movement <- read.csv(paste0(data_path,"movement_data_HLM_-14toEnd_allparticipants.csv"))
names(df_hlm_movement)[6:13] <- paste0("movement_",names(df_hlm_movement)[6:13])

df_hlm_time <- read.csv(paste0(data_path,"time_data_HLM_-14toEnd_allparticipants.csv"))
names(df_hlm_time)[6:13] <- paste0("time_",names(df_hlm_time)[6:13])

df_hlm_mood <- read.csv(paste0(data_path,"mood_data_HLM_-14toEnd_allparticipants.csv"))

df_hlm <- Reduce(merge, list(df_hlm_movement,df_hlm_time,df_hlm_mood))

df_demo <- read.csv(paste0(data_path,"/demo_info.csv"))%>%
  mutate(Allocation = if_else(Allocation == "B","Lithium", "Control"))

participant_num <- length(unique(df_hlm$ID))

df_hlm$Allocation[df_hlm$Allocation == "A"] <- "control"
df_hlm$Allocation[df_hlm$Allocation == "B"] <- "lithium"

# mood parameters from 0~4 to 1~5
df_hlm$pos_mean <- df_hlm$pos_mean+5
df_hlm$neg_mean <- df_hlm$neg_mean+5

# mood parameters need log transformation
df_hlm$posLog <- log(df_hlm$pos_mean+1)
df_hlm$negLog <- log(df_hlm$neg_mean+1)

#compute relative amplitude
df_hlm <- df_hlm%>%
  rowwise()%>%
  mutate(rel_amplitude  = (movement_M10_mean - movement_L5_mean)/(movement_M10_mean + movement_L5_mean))

#add season into df
df_hlm <- merge(df_hlm,df_demo[,c("ID","Season","Diagnosis")], by = "ID")

data availability

group_N <- count(unique(df_hlm%>%
  select(ID,Allocation)),Allocation)

# data availability by week (rolling mean)
ymax = max(group_N$n)+4
p1 <- ggplot(df_hlm%>%
  filter(!is.na(movement_M10_mean))%>%
  group_by(Allocation)%>%
  count(day),
  aes(x = day, y = n, color = Allocation))+
  geom_line(size=1.5)+
  geom_segment(x = -14, xend = -14, y = ymax, yend = ymax-4, color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  geom_segment(x = 28, xend = 28, y = ymax, yend = ymax-4, color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  labs(title = "Movement data availability",
       y = "available N")+
  ylim(0,max(group_N$n)+1)+
  scale_color_manual(values = orangeNavyPalette)+
  scale_x_continuous(breaks = seq(-14,56,by = 7))+
  theme(legend.position = "bottom",
        panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.

p2 <- ggplot(df_hlm%>%
  filter(!is.na(pos_mean))%>%
  group_by(Allocation)%>%
  count(day),
  aes(x = day, y = n, color = Allocation))+
  geom_line(size=1.5)+
  geom_segment(x = -14, xend = -14, y = ymax, yend = ymax-4, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  geom_segment(x = 28, xend = 28, y = ymax, yend = ymax-4, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  labs(title = "Mood data availability",
       y = "available N")+
  ylim(0,max(group_N$n)+1)+
  scale_color_manual(values = orangeNavyPalette)+
  scale_x_continuous(breaks = seq(-14,56,by = 7))+
  theme(legend.position = "none",
        panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())

p3 <- ggplot( df_hlm%>%
  filter(!is.na(movement_M10_vmu))%>%
  group_by(Allocation)%>%
  count(day),
  aes(x = day, y = n, color = Allocation))+
  geom_line(size=1.5)+
  geom_segment(x = -14, xend = -14, y = ymax, yend = ymax-4, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  geom_segment(x = 28, xend = 28, y = ymax, yend = ymax-4, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  labs(title = "Movement model output availability",
       y = "available N")+
  ylim(0,max(group_N$n)+1)+
  scale_color_manual(values = orangeNavyPalette)+
  scale_x_continuous(breaks = seq(-14,56,by = 7))+
  theme(legend.position = "none",
        panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())

p4 <- ggplot(df_hlm%>%
  filter(!is.na(pos_vmu))%>%
  group_by(Allocation)%>%
  count(day), aes(x = day, y = n, color = Allocation))+
  geom_line(size=1.5)+
  geom_segment(x = -14, xend = -14, y = ymax, yend = ymax -4, color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  geom_segment(x = 28, xend = 28, y = ymax, yend = ymax -4, color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  labs(title = "Mood model output availability",
       y = "available N")+
  ylim(0,max(group_N$n)+1)+
  scale_color_manual(values = orangeNavyPalette)+
  scale_x_continuous(breaks = seq(-14,56,by = 7))+
  theme(legend.position = "none",
        panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())

group_N <- count(unique(df_hlm%>%
  select(ID,Allocation)),Allocation)

# data availability by week (rolling mean)
ymax = max(group_N$n)+4
p1a <- ggplot(df_hlm%>%
  filter(!is.na(movement_M10_mean))%>%
  group_by(Allocation)%>%
  count(day),
  aes(x = day, y = n, color = Allocation))+
  geom_line(size=1.5)+
  geom_segment(x = -14, xend = -14, y = ymax, yend = ymax-6, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  geom_segment(x = 28, xend = 28, y = ymax, yend = ymax-6, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  labs(title = "Movement data availability",
       y = "available N")+
  ylim(0,max(group_N$n)+1)+
  scale_color_manual(values = orangeNavyPalette)+
  scale_x_continuous(breaks = seq(-14,56,by = 7))+
  theme(legend.position = "bottom",
        panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())

p1b <- ggplot(df_hlm%>%
  filter(!is.na(pos_mean))%>%
  group_by(Allocation)%>%
  count(day),
  aes(x = day, y = n, color = Allocation))+
  geom_line(size=1.5)+
  geom_segment(x = -14, xend = -14, y = ymax, yend = ymax-6, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  geom_segment(x = 28, xend = 28, y = ymax, yend = ymax-6, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  labs(title = "Mood data availability",
       y = "available N")+
  ylim(0,max(group_N$n)+1)+
  scale_color_manual(values = orangeNavyPalette)+
  scale_x_continuous(breaks = seq(-14,56,by = 7))+
  theme(legend.position = "none",
        panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())

p1c <- ggplot( df_hlm%>%
  filter(!is.na(movement_M10_vmu))%>%
  group_by(Allocation)%>%
  count(day),
  aes(x = day, y = n, color = Allocation))+
  geom_line(size=1.5)+
  geom_segment(x = -14, xend = -14, y = ymax, yend = ymax-6, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  geom_segment(x = 28, xend = 28, y = ymax, yend = ymax-6, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  labs(title = "Movement model output availability",
       y = "available N")+
  ylim(0,max(group_N$n)+1)+
  scale_color_manual(values = orangeNavyPalette)+
  scale_x_continuous(breaks = seq(-14,56,by = 7))+
  theme(legend.position = "none",
        panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())

p1d <- ggplot(df_hlm%>%
  filter(!is.na(pos_vmu))%>%
  group_by(Allocation)%>%
  count(day), aes(x = day, y = n, color = Allocation))+
  geom_line(size=1.5)+
  geom_segment(x = -14, xend = -14, y = ymax, yend = ymax -6, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  geom_segment(x = 28, xend = 28, y = ymax, yend = ymax -6, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  labs(title = "Mood model output availability",
       y = "available N")+
  ylim(0,max(group_N$n)+1)+
  scale_color_manual(values = orangeNavyPalette)+
  scale_x_continuous(breaks = seq(-14,56,by = 7))+
  theme(legend.position = "none",
        panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())

(p1a+p1b)/(p1c+p1d)+ 
  plot_annotation(tag_levels = 'A')

Data availability; grey arrows represent the period of time considered by this study.

  ggsave("../figs/S1_data_availability.png",width = 8 ,height=6,dpi = 300)

df_hlm <- df_hlm%>%filter(day >= -14 & day <= 28)

density plot of movement parameters

p1d <- ggdensity(df_hlm$movement_M10_mean)+labs(title = "M10 mean",x="")+
  scale_x_continuous(breaks = seq(25,200,by = 75))
p2d <- ggdensity(df_hlm$movement_M10_s)+labs(title = "M10 noise",x="")
p3d <- ggdensity(df_hlm$movement_M10_vmu)+labs(title = "M10 volatility",x="")
p4d <- ggdensity(df_hlm$movement_M10_mu)+labs(title = "M10 mu",x="")
p5d <- ggdensity(df_hlm$movement_L5_mean)+labs(title = "L5 mean",x="")
p6d <- ggdensity(df_hlm$movement_L5_s)+labs(title = "L5 noise",x="")
p7d <- ggdensity(df_hlm$movement_L5_vmu)+labs(title = "L5 volatility",x="")
p8d <- ggdensity(df_hlm$movement_L5_mu)+labs(title = "L5 mu",x="")+
  scale_x_continuous(breaks = seq(0.1,0.5,by = 0.1))

(p1d|p2d|p3d|p4d)/(p5d|p6d|p7d|p8d)+
  plot_annotation(title = "distribution of movement parameters")

## Warning: Removed 695 rows containing non-finite values (`stat_density()`).

## Warning: Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).

## Warning: Removed 695 rows containing non-finite values (`stat_density()`).

## Warning: Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).

 ggsave("../figs/params_distribution_movement.png",width = 8 ,height=5,dpi = 300)

## Warning: Removed 695 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).

## Warning: Removed 695 rows containing non-finite values (`stat_density()`).

## Warning: Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).

density plot of movement onset time parameters

p1d <- ggdensity(df_hlm$time_M10_mean, linetype = "dashed")+labs(title = "M10 time mean",x="")
p2d <- ggdensity(df_hlm$time_M10_s, linetype = "dashed")+labs(title = "M10 time noise",x="")
p3d <- ggdensity(df_hlm$time_M10_vmu, linetype = "dashed")+labs(title = "M10 time volatility",x="")
p4d <- ggdensity(df_hlm$time_M10_mu, linetype = "dashed")+labs(title = "M10 time mu",x="")
p5d <- ggdensity(df_hlm$time_L5_mean, linetype = "dashed")+labs(title = "L5 time mean",x="")
p6d <- ggdensity(df_hlm$time_L5_s, linetype = "dashed")+labs(title = "L5 time noise",x="")
p7d <- ggdensity(df_hlm$time_L5_vmu, linetype = "dashed")+labs(title = "L5 time volatility",x="")
p8d <- ggdensity(df_hlm$time_L5_mu, linetype = "dashed")+labs(title = "L5 time mu",x="")+
  scale_x_continuous(breaks = seq(0,1,by = 0.3))

(p1d|p2d|p3d|p4d)/(p5d|p6d|p7d|p8d)+
  plot_annotation(title = "distribution of onset time parameters")

## Warning: Removed 695 rows containing non-finite values (`stat_density()`).

## Warning: Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).

## Warning: Removed 695 rows containing non-finite values (`stat_density()`).

## Warning: Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).

 ggsave("../figs/params_distribution_onsetTime.png",width = 8 ,height=5,dpi = 300)

## Warning: Removed 695 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).

## Warning: Removed 695 rows containing non-finite values (`stat_density()`).

## Warning: Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).

density plot of mood

p1d <- ggdensity(df_hlm$posLog)+labs(title = "log PA")
p2d <- ggdensity(df_hlm$negLog)+labs(title = "log NA")
p3d <- ggdensity(df_hlm$pos_vmu)+labs(title = "PA vmu")
p4d <- ggdensity(df_hlm$neg_vmu)+labs(title = "NA vmu")

(p1d|p2d)/(p3d|p4d)

## Warning: Removed 577 rows containing non-finite values (`stat_density()`).
## Removed 577 rows containing non-finite values (`stat_density()`).

## Warning: Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).

df_hlm$phase[df_hlm$day < 1] <- "pre"
df_hlm$phase[df_hlm$day >= 1 & df_hlm$day <= 7] <- "post1"
df_hlm$phase[df_hlm$day >= 8 & df_hlm$day <= 14] <- "post2"
df_hlm$phase[df_hlm$day >= 15 & df_hlm$day <= 21] <- "post3"
df_hlm$phase[df_hlm$day >= 22 & df_hlm$day <= 28] <- "post4"

df_hlm$phase <- factor(df_hlm$phase, levels = c("pre","post1","post2","post3","post4")) #,"post5","post6"

df_hlm <- df_hlm%>%
  relocate(phase, .after= day)

correlation among variables

#correlation heatmap

find.corr <- function(df,a,b){
  #df is the dataframe
  #a(IV) and b(DV) are the column index to which the model is fitted
  iv <- names(df)[a]
  dv <- names(df)[b]
  my.formula = as.formula(paste0("scale(",iv,") ~ scale(",dv,") + (1 + scale(",dv,") |ID)"))
  my.lme <- summary(lmer(my.formula,df))
  return(c(my.lme$coefficients[2,1],my.lme$coefficients[2,5])) #beta, p
}

df_hlm_corr2 <- df_hlm%>%
  select(ID,day, Allocation,movement_M10_mean,
         time_M10_mean,posLog,negLog,
         movement_M10_vmu,
         time_M10_vmu,
         pos_vmu,neg_vmu,
         movement_M10_s,
         time_M10_s,
         pos_s,neg_s)%>%
  rename(M10 = movement_M10_mean,
         timeM10 = time_M10_mean,
         logPA = posLog,
         logNA = negLog,
         timeM10_vmu = time_M10_vmu,
         M10_vmu = movement_M10_vmu,
         PA_vmu = pos_vmu,
         NA_vmu = neg_vmu,
         M10_s = movement_M10_s,
         timeM10_s = time_M10_s,
         PA_s = pos_s,
         NA_s = neg_s)

num.row2 = ncol(df_hlm_corr2)-3

my.corr.mat2 <- matrix(nrow = num.row2, ncol = num.row2)
my.corr.pmat2 <- matrix(nrow = num.row2, ncol = num.row2)

rownames(my.corr.mat2) <- names(df_hlm_corr2)[4:(num.row2+3)]
colnames(my.corr.mat2) <- names(df_hlm_corr2)[4:(num.row2+3)]
rownames(my.corr.pmat2) <- names(df_hlm_corr2)[4:(num.row2+3)]
colnames(my.corr.pmat2) <- names(df_hlm_corr2)[4:(num.row2+3)]

for (i in 1:num.row2){
  my.corr.mat2[i,i] = 1
  my.corr.pmat2[i,i] = 0
}

for (a in 4:(num.row2+3)){
  for (b in 4:(num.row2+3)){
    if (a != b){
      temp_out <- find.corr(df_hlm_corr2,a,b)
      my.corr.mat2[a-3,b-3] = temp_out[1]
      my.corr.pmat2[a-3,b-3] = temp_out[2]
    }
  }
}

## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00201368 (tol = 0.002, component 1)

## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')

ggcorrplot::ggcorrplot(my.corr.mat2,p.mat = my.corr.pmat2,insig = "blank",colors=c("orange", "white", "navy"),lab = TRUE, digits = 1)

round(my.corr.mat2,3)

##                M10 timeM10  logPA  logNA M10_vmu timeM10_vmu PA_vmu NA_vmu
## M10          1.000  -0.031  0.194 -0.076  -0.002       0.001 -0.066 -0.029
## timeM10     -0.045   1.000 -0.061  0.047   0.059       0.071  0.003  0.041
## logPA        0.150  -0.007  1.000 -0.172   0.025       0.080  0.064  0.003
## logNA       -0.061   0.019 -0.188  1.000   0.153       0.156 -0.064  0.037
## M10_vmu     -0.007   0.037  0.058  0.178   1.000       0.592  0.283  0.331
## timeM10_vmu  0.005   0.054  0.122  0.153   0.593       1.000  0.273  0.306
## PA_vmu      -0.055   0.015  0.035 -0.054   0.143       0.036  1.000  0.245
## NA_vmu      -0.018   0.057 -0.049  0.075   0.305       0.088  0.401  1.000
## M10_s        0.022   0.038  0.040 -0.024   0.313       0.196 -0.007 -0.008
## timeM10_s   -0.030  -0.048 -0.035  0.053  -0.144      -0.162  0.045  0.005
## PA_s        -0.028   0.049  0.044  0.105   0.130       0.294 -0.053  0.077
## NA_s         0.014   0.041  0.009  0.091   0.148       0.245  0.191  0.260
##              M10_s timeM10_s   PA_s  NA_s
## M10         -0.005    -0.024  0.001 0.005
## timeM10      0.066    -0.057  0.073 0.065
## logPA        0.072    -0.060 -0.010 0.000
## logNA        0.010     0.029  0.040 0.053
## M10_vmu      0.390    -0.199  0.114 0.248
## timeM10_vmu  0.286    -0.204  0.204 0.350
## PA_vmu      -0.033    -0.035 -0.178 0.097
## NA_vmu       0.023     0.049 -0.096 0.256
## M10_s        1.000    -0.031  0.178 0.053
## timeM10_s   -0.052     1.000 -0.030 0.002
## PA_s         0.133    -0.020  1.000 0.332
## NA_s         0.083    -0.073  0.203 1.000

round(my.corr.pmat2,3)

##               M10 timeM10 logPA logNA M10_vmu timeM10_vmu PA_vmu NA_vmu M10_s
## M10         0.000   0.443 0.007 0.167   0.960       0.987  0.122  0.423 0.903
## timeM10     0.453   0.000 0.274 0.445   0.141       0.077  0.947  0.309 0.122
## logPA       0.003   0.870 0.000 0.022   0.683       0.212  0.398  0.954 0.251
## logNA       0.114   0.579 0.005 0.000   0.007       0.019  0.255  0.330 0.814
## M10_vmu     0.857   0.269 0.390 0.004   0.000       0.000  0.014  0.002 0.000
## timeM10_vmu 0.897   0.102 0.055 0.028   0.000       0.000  0.035  0.004 0.002
## PA_vmu      0.077   0.600 0.564 0.298   0.069       0.711  0.000  0.007 0.673
## NA_vmu      0.614   0.087 0.514 0.325   0.002       0.400  0.001  0.000 0.771
## M10_s       0.502   0.229 0.487 0.577   0.000       0.005  0.924  0.898 0.000
## timeM10_s   0.411   0.256 0.540 0.408   0.065       0.098  0.638  0.942 0.462
## PA_s        0.609   0.097 0.378 0.180   0.058       0.002  0.621  0.314 0.040
## NA_s        0.727   0.137 0.828 0.041   0.040       0.003  0.034  0.003 0.296
##             timeM10_s  PA_s  NA_s
## M10             0.512 0.983 0.930
## timeM10         0.242 0.140 0.209
## logPA           0.309 0.860 0.999
## logNA           0.621 0.437 0.218
## M10_vmu         0.079 0.351 0.021
## timeM10_vmu     0.095 0.158 0.002
## PA_vmu          0.682 0.102 0.296
## NA_vmu          0.577 0.384 0.032
## M10_s           0.739 0.025 0.530
## timeM10_s       0.000 0.719 0.980
## PA_s            0.791 0.000 0.000
## NA_s            0.303 0.010 0.000

ggsave("../figs/S2_corr_plot.png",width = 6 ,height=6,dpi = 300)

pre-randomisation check MANOVA

# check M10
df1 <- df_hlm%>%
  dplyr::filter(phase == "pre")%>%
  group_by(Allocation,ID)%>%
  summarise_at(vars(movement_L5_mean:time_M10_mu), ~mean(.x, na.rm = TRUE))

m01 <- manova(cbind(movement_M10_mean,movement_M10_vmu,time_M10_mean,time_M10_vmu) ~ Allocation, data =df1)
summary.aov(m01)

##  Response movement_M10_mean :
##             Df  Sum Sq Mean Sq F value Pr(>F)
## Allocation   1    16.9   16.90  0.0476 0.8289
## Residuals   29 10300.2  355.18               
## 
##  Response movement_M10_vmu :
##             Df Sum Sq Mean Sq F value Pr(>F)
## Allocation   1 0.4642 0.46422  2.6364 0.1153
## Residuals   29 5.1064 0.17608               
## 
##  Response time_M10_mean :
##             Df  Sum Sq Mean Sq F value Pr(>F)
## Allocation   1   0.922  0.9217  0.1868 0.6688
## Residuals   29 143.100  4.9345               
## 
##  Response time_M10_vmu :
##             Df  Sum Sq  Mean Sq F value Pr(>F)
## Allocation   1 0.00436 0.004359  0.0484 0.8273
## Residuals   29 2.60987 0.089996               
## 
## 3 observations deleted due to missingness

None of the parameters of interest were different pre-randomization.

#linear mixed model ## mood ### raw value

model_formula_mean <- posLog ~ Allocation * phase + Age + Gender + (1|ID)

df_hlm_lmer_mood_mean <- df_hlm%>%
  select(ID,Allocation,phase,day,posLog,Age,Gender,Season)

mod_mood_mean_pos <- lmer(model_formula_mean,data = df_hlm_lmer_mood_mean)
kable(round(summary(mod_mood_mean_pos)$coefficients,3))

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	2.439	0.165	30.858	14.787	0.000
Allocationlithium	0.075	0.102	36.666	0.736	0.466
phasepost1	-0.114	0.039	843.496	-2.930	0.003
phasepost2	-0.028	0.040	845.510	-0.713	0.476
phasepost3	0.021	0.040	845.093	0.519	0.604
phasepost4	0.033	0.040	845.056	0.824	0.410
Age	0.002	0.004	29.512	0.454	0.653
GenderM	-0.277	0.095	29.882	-2.904	0.007
Allocationlithium:phasepost1	0.019	0.053	844.156	0.358	0.721
Allocationlithium:phasepost2	-0.076	0.055	844.619	-1.384	0.167
Allocationlithium:phasepost3	-0.092	0.056	844.998	-1.647	0.100
Allocationlithium:phasepost4	-0.144	0.057	845.539	-2.534	0.011

lme.dscore(mod_mood_mean_pos,df_hlm_lmer_mood_mean,type = "lme4")%>%
  as.data.frame()%>%
  kable()

	t	df	d
Allocationlithium	0.7361081	36.66568	0.2431318
phasepost1	-2.9295452	843.49614	-0.2017384
phasepost2	-0.7132416	845.50969	-0.0490577
phasepost3	0.5190409	845.09262	0.0357091
phasepost4	0.8240219	845.05637	0.0566925
Age	0.4540456	29.51158	0.1671603
GenderM	-2.9044919	29.88190	-1.0626643
Allocationlithium:phasepost1	0.3577954	844.15572	0.0246294
Allocationlithium:phasepost2	-1.3839384	844.61873	-0.0952393
Allocationlithium:phasepost3	-1.6466073	844.99783	-0.1132902
Allocationlithium:phasepost4	-2.5339154	845.53869	-0.1742831

ee1 <- as.data.frame(Effect(c("phase", "Allocation"),mod_mood_mean_pos))

stat.test <- as.data.frame(round(summary(mod_mood_mean_pos)$coefficients,3))[9:12,]%>%
  rename(p_val = `Pr(>|t|)`)

stat.test$group1 <- c("post1","post2","post3","post4")
stat.test$group2 <- c("post1","post2","post3","post4")

stat.test <- left_join(stat.test,ee1[ee1$phase != "pre" & ee1$Allocation == "control",]%>%
  rename(group1 = phase))%>%
  mutate(y.position = ee1$upper[1])%>%
  mutate(p_sig = if_else(p_val < 0.001,"***", if_else(p_val < 0.01, "**", if_else(p_val < 0.05, "*", "ns")) ))

## Joining, by = "group1"

p8 <- ggplot(ee1, aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank(),
        axis.text.x = element_text(angle=45))+
  labs(title = "",y="mood level")+ 
  stat_pvalue_manual(stat.test, label = "p_sig",vjust = -1)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))
p8

ggsave("../figs/hlm_mood_val.png",width = 4 ,height=3,dpi = 300)

variability

model_formula <- val ~ Allocation * var_type * phase + Age + Gender + (var_type|ID)

df_hlm_lmer_mood_variability <- df_hlm%>%
  select(ID,Allocation,day,phase,pos_vmu,pos_s,neg_s,Age,Gender,posLog)%>%
  pivot_longer(c(pos_vmu,pos_s),names_pattern =  "(.*)_(.*)",names_to = c("mood_type","var_type"),values_to = "val")%>%
  mutate(var_type = if_else(var_type == "vmu","volatility","noise"))

df_hlm_lmer_mood_variability$var_type <- factor(df_hlm_lmer_mood_variability$var_type)

df_hlm_lmer_mood_variability$var_type <- relevel(df_hlm_lmer_mood_variability$var_type, ref = "volatility")

mod_pos_vmu <- lmer(model_formula,data = df_hlm_lmer_mood_variability)
kable(round(summary(mod_pos_vmu)$coefficients,3))

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	-3.658	0.300	54.622	-12.198	0.000
Allocationlithium	-0.078	0.321	34.460	-0.244	0.809
var_typenoise	2.688	0.268	36.165	10.012	0.000
phasepost1	-0.522	0.079	2402.420	-6.614	0.000
phasepost2	-0.747	0.080	2402.793	-9.341	0.000
phasepost3	-0.706	0.081	2403.788	-8.721	0.000
phasepost4	-0.618	0.081	2403.788	-7.635	0.000
Age	-0.003	0.005	29.865	-0.597	0.555
GenderM	-0.310	0.117	30.016	-2.652	0.013
Allocationlithium:var_typenoise	0.052	0.359	36.037	0.144	0.886
Allocationlithium:phasepost1	0.062	0.105	2402.447	0.586	0.558
Allocationlithium:phasepost2	0.227	0.106	2402.661	2.145	0.032
Allocationlithium:phasepost3	0.554	0.107	2403.242	5.194	0.000
Allocationlithium:phasepost4	0.742	0.107	2403.242	6.964	0.000
var_typenoise:phasepost1	0.463	0.112	2403.171	4.144	0.000
var_typenoise:phasepost2	0.611	0.113	2404.255	5.403	0.000
var_typenoise:phasepost3	0.607	0.114	2407.094	5.298	0.000
var_typenoise:phasepost4	0.475	0.114	2407.094	4.144	0.000
Allocationlithium:var_typenoise:phasepost1	-0.092	0.149	2403.255	-0.616	0.538
Allocationlithium:var_typenoise:phasepost2	-0.227	0.150	2403.887	-1.515	0.130
Allocationlithium:var_typenoise:phasepost3	-0.648	0.151	2405.567	-4.301	0.000
Allocationlithium:var_typenoise:phasepost4	-0.793	0.151	2405.567	-5.259	0.000

ee2 <- as.data.frame(Effect(c("phase", "Allocation","var_type"),mod_pos_vmu))%>%
  filter(var_type == "volatility")

stat.test_pos_vmu <- as.data.frame(round(summary(mod_pos_vmu)$coefficients,3))[11:14,]

stat.test <- stat.test_pos_vmu%>%rename(p_val = 'Pr(>|t|)')
stat.test$group1 <- c("post1","post2","post3","post4")
stat.test$group2 <- c("post1","post2","post3","post4")

stat.test <- left_join(stat.test,ee2[ee2$phase != "pre" & ee2$Allocation == "lithium",]%>%
  rename(group1 = phase), by = c("group1"))%>%
  mutate(y.position = ee2$upper[1])%>%
  mutate(p_sig = if_else(p_val < 0.001,"***", if_else(p_val < 0.01, "**", if_else(p_val < 0.05, "*", "ns")) ))


#plot volatility first
p9 <- ggplot(ee2%>%filter(var_type == "volatility"), aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank())+
  labs(title = "",y="mood volatility")+
  scale_y_continuous(limits = c(-5.4,-3),breaks = seq(-5.4,-3, by = 0.4))+
  stat_pvalue_manual(stat.test%>%filter(var_type == "volatility"),label = "p_sig",vjust = -1)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))
p9

ggsave("../figs/mood_variability.png",width = 4 ,height=3,dpi = 300)

movement

relative amplitude

model_formula_amp <- rel_amplitude ~ Allocation * phase + posLog + pos_vmu + Age + Gender + (1|ID)
mod_relative_amplitude <- lmer(model_formula_amp,data = df_hlm)
kable(round(summary(mod_relative_amplitude)$coefficients,3))

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	0.770	0.035	110.921	22.316	0.000
Allocationlithium	-0.018	0.016	49.303	-1.129	0.264
phasepost1	0.005	0.009	538.256	0.488	0.626
phasepost2	0.014	0.010	543.395	1.367	0.172
phasepost3	0.012	0.009	543.427	1.342	0.180
phasepost4	0.018	0.010	550.841	1.703	0.089
posLog	0.037	0.009	533.222	4.341	0.000
pos_vmu	-0.003	0.003	536.463	-1.281	0.201
Age	0.000	0.001	30.355	-0.851	0.402
GenderM	0.004	0.014	32.775	0.260	0.797
Allocationlithium:phasepost1	-0.008	0.014	539.910	-0.557	0.578
Allocationlithium:phasepost2	0.000	0.014	541.911	0.013	0.990
Allocationlithium:phasepost3	-0.006	0.013	545.726	-0.411	0.681
Allocationlithium:phasepost4	-0.017	0.015	553.843	-1.128	0.260

plot_model(mod_relative_amplitude, type = "pred",terms = c("phase", "Allocation"))

raw value

model_formula_mean <- val ~ Allocation * movement_type * phase + posLog + pos_vmu  + Age + Gender + (movement_type|ID)

df_hlm_lmer_movement_mean <- df_hlm%>%
  select(ID,Allocation,phase,day,movement_L5_mean,movement_M10_mean,Age,Gender,posLog, negLog, pos_vmu,Season)%>%
  pivot_longer(c(movement_L5_mean,movement_M10_mean),names_pattern =  "movement_(.*)_mean",names_to = c("movement_type"),values_to = "val")

df_hlm_lmer_movement_mean$movement_type <- factor(df_hlm_lmer_movement_mean$movement_type)
df_hlm_lmer_movement_mean$movement_type <- relevel(df_hlm_lmer_movement_mean$movement_type, ref = "M10")

mod_movement_mean_M10 <- lmer(model_formula_mean,data = df_hlm_lmer_movement_mean)

## boundary (singular) fit: see help('isSingular')

anova(mod_movement_mean_M10)

## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq Mean Sq NumDF   DenDF  F value    Pr(>F)
## Allocation                        663     663     1   31.82   2.9724  0.094400
## movement_type                   82158   82158     1   32.16 368.0827 < 2.2e-16
## phase                             332      83     4 1091.39   0.3721  0.828647
## posLog                           2120    2120     1  757.53   9.4960  0.002134
## pos_vmu                           408     408     1  804.91   1.8274  0.176816
## Age                                22      22     1  237.17   0.0972  0.755540
## Gender                            302     302     1  266.39   1.3541  0.245596
## Allocation:movement_type          595     595     1   32.13   2.6644  0.112381
## Allocation:phase                 1511     378     4 1093.41   1.6927  0.149344
## movement_type:phase               445     111     4 1094.08   0.4988  0.736646
## Allocation:movement_type:phase   1974     494     4 1096.48   2.2111  0.065843
##                                   
## Allocation                     .  
## movement_type                  ***
## phase                             
## posLog                         ** 
## pos_vmu                           
## Age                               
## Gender                            
## Allocation:movement_type          
## Allocation:phase                  
## movement_type:phase               
## Allocation:movement_type:phase .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

kable(round(summary(mod_movement_mean_M10)$coefficients,3))

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	44.642	6.399	171.727	6.977	0.000
Allocationlithium	-0.879	5.809	42.440	-0.151	0.880
movement_typeL5	-53.929	4.636	48.866	-11.632	0.000
phasepost1	2.363	2.643	1097.183	0.894	0.372
phasepost2	5.807	2.867	1105.181	2.026	0.043
phasepost3	5.127	2.652	1108.042	1.933	0.053
phasepost4	6.989	2.963	1115.918	2.359	0.019
posLog	4.494	1.458	757.526	3.082	0.002
pos_vmu	-0.613	0.454	804.907	-1.352	0.177
Age	-0.017	0.054	237.166	-0.312	0.756
GenderM	1.564	1.344	266.386	1.164	0.246
Allocationlithium:movement_typeL5	0.366	6.329	52.646	0.058	0.954
Allocationlithium:phasepost1	-8.375	3.917	1103.919	-2.138	0.033
Allocationlithium:phasepost2	-8.055	3.973	1106.937	-2.028	0.043
Allocationlithium:phasepost3	-12.124	3.849	1111.723	-3.150	0.002
Allocationlithium:phasepost4	-13.768	4.187	1119.807	-3.288	0.001
movement_typeL5:phasepost1	-1.522	3.704	1095.622	-0.411	0.681
movement_typeL5:phasepost2	-5.849	3.970	1095.939	-1.473	0.141
movement_typeL5:phasepost3	-5.839	3.670	1096.865	-1.591	0.112
movement_typeL5:phasepost4	-7.781	4.053	1089.207	-1.920	0.055
Allocationlithium:movement_typeL5:phasepost1	7.817	5.474	1097.558	1.428	0.154
Allocationlithium:movement_typeL5:phasepost2	8.191	5.525	1097.678	1.482	0.139
Allocationlithium:movement_typeL5:phasepost3	12.797	5.315	1097.084	2.408	0.016
Allocationlithium:movement_typeL5:phasepost4	14.929	5.720	1091.567	2.610	0.009

df_hlm_lmer_movement_mean$movement_type <- factor(df_hlm_lmer_movement_mean$movement_type)
df_hlm_lmer_movement_mean$movement_type <- relevel(df_hlm_lmer_movement_mean$movement_type, ref = "L5")

mod_movement_mean_L5 <- lmer(model_formula_mean,data = df_hlm_lmer_movement_mean)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00627183 (tol = 0.002, component 1)

kable(round(summary(mod_movement_mean_L5)$coefficients,3))

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	-9.286	5.041	757.080	-1.842	0.066
Allocationlithium	-0.513	2.474	1095.870	-0.207	0.836
movement_typeM10	53.929	4.637	48.839	11.630	0.000
phasepost1	0.840	2.625	1095.144	0.320	0.749
phasepost2	-0.042	2.780	1095.401	-0.015	0.988
phasepost3	-0.712	2.561	1095.736	-0.278	0.781
phasepost4	-0.793	2.795	1095.070	-0.284	0.777
posLog	4.494	1.458	757.535	3.082	0.002
pos_vmu	-0.613	0.454	804.924	-1.352	0.177
Age	-0.017	0.054	237.231	-0.313	0.755
GenderM	1.565	1.344	266.457	1.165	0.245
Allocationlithium:movement_typeM10	-0.366	6.330	52.616	-0.058	0.954
Allocationlithium:phasepost1	-0.557	3.825	1095.232	-0.146	0.884
Allocationlithium:phasepost2	0.135	3.848	1095.589	0.035	0.972
Allocationlithium:phasepost3	0.673	3.691	1096.204	0.182	0.855
Allocationlithium:phasepost4	1.161	3.941	1094.894	0.295	0.768
movement_typeM10:phasepost1	1.523	3.704	1095.618	0.411	0.681
movement_typeM10:phasepost2	5.849	3.970	1095.922	1.473	0.141
movement_typeM10:phasepost3	5.839	3.670	1096.843	1.591	0.112
movement_typeM10:phasepost4	7.781	4.053	1089.167	1.920	0.055
Allocationlithium:movement_typeM10:phasepost1	-7.818	5.474	1097.546	-1.428	0.154
Allocationlithium:movement_typeM10:phasepost2	-8.191	5.525	1097.659	-1.482	0.139
Allocationlithium:movement_typeM10:phasepost3	-12.797	5.315	1097.056	-2.408	0.016
Allocationlithium:movement_typeM10:phasepost4	-14.929	5.720	1091.525	-2.610	0.009

plot_model(mod_movement_mean_M10, type = "pred", terms = c("phase", "Allocation","movement_type"))+
  scale_color_manual(values=orangeNavyPalette)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.5),
         panel.grid.minor.y = element_blank())+
  labs(title = "predicted value of movement level",y="movement level")+
  scale_y_continuous(breaks = seq(-10, 70, by = 10))

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

df_hlm_lmer_movement_mean$Allocation <- factor(df_hlm_lmer_movement_mean$Allocation, levels = c("lithium","control"))
kable(round(summary(lmer(model_formula_mean,data = df_hlm_lmer_movement_mean))$coefficients,3))

## Warning: Model failed to converge with 1 negative eigenvalue: -3.9e-02

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	-9.800	5.068	793.682	-1.934	0.053
Allocationcontrol	0.513	2.474	1096.082	0.207	0.836
movement_typeM10	53.563	4.308	57.573	12.434	0.000
phasepost1	0.282	2.798	1095.388	0.101	0.920
phasepost2	0.094	2.675	1095.647	0.035	0.972
phasepost3	-0.039	2.641	1095.405	-0.015	0.988
phasepost4	0.369	2.773	1094.914	0.133	0.894
posLog	4.494	1.458	757.525	3.082	0.002
pos_vmu	-0.613	0.454	804.906	-1.352	0.177
Age	-0.017	0.054	237.165	-0.312	0.756
GenderM	1.564	1.344	266.385	1.164	0.246
Allocationcontrol:movement_typeM10	0.366	6.329	52.644	0.058	0.954
Allocationcontrol:phasepost1	0.558	3.825	1095.233	0.146	0.884
Allocationcontrol:phasepost2	-0.135	3.848	1095.589	-0.035	0.972
Allocationcontrol:phasepost3	-0.673	3.691	1096.205	-0.182	0.855
Allocationcontrol:phasepost4	-1.161	3.941	1094.893	-0.295	0.768
movement_typeM10:phasepost1	-6.295	4.027	1095.519	-1.563	0.118
movement_typeM10:phasepost2	-2.342	3.835	1096.171	-0.611	0.541
movement_typeM10:phasepost3	-6.958	3.835	1092.424	-1.814	0.070
movement_typeM10:phasepost4	-7.148	4.039	1094.587	-1.770	0.077
Allocationcontrol:movement_typeM10:phasepost1	7.817	5.474	1097.559	1.428	0.154
Allocationcontrol:movement_typeM10:phasepost2	8.191	5.525	1097.678	1.482	0.139
Allocationcontrol:movement_typeM10:phasepost3	12.797	5.315	1097.084	2.408	0.016
Allocationcontrol:movement_typeM10:phasepost4	14.929	5.720	1091.567	2.610	0.009

df_hlm_lmer_movement_mean$Allocation <- relevel(df_hlm_lmer_movement_mean$Allocation, ref = "control")

relevel to different weeks to show group contrast

df_hlm_lmer_movement_mean$movement_type <- relevel(df_hlm_lmer_movement_mean$movement_type, ref = "M10")
df_hlm_lmer_movement_mean$phase <- relevel(df_hlm_lmer_movement_mean$phase, ref = "post1")
kable(round(summary(lmer(model_formula_mean,data = df_hlm_lmer_movement_mean))$coefficients,3))

## boundary (singular) fit: see help('isSingular')

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	47.004	6.440	174.854	7.298	0.000
Allocationlithium	-9.254	6.089	51.066	-1.520	0.135
movement_typeL5	-55.451	4.848	58.069	-11.437	0.000
phasepre	-2.363	2.643	1097.183	-0.894	0.372
phasepost2	3.444	3.024	1105.717	1.139	0.255
phasepost3	2.765	2.851	1108.841	0.970	0.332
phasepost4	4.626	3.138	1117.351	1.474	0.141
posLog	4.494	1.458	757.526	3.082	0.002
pos_vmu	-0.613	0.454	804.907	-1.352	0.177
Age	-0.017	0.054	237.166	-0.312	0.756
GenderM	1.564	1.344	266.386	1.164	0.246
Allocationlithium:movement_typeL5	8.183	6.800	69.303	1.203	0.233
Allocationlithium:phasepre	8.375	3.917	1103.919	2.138	0.033
Allocationlithium:phasepost2	0.320	4.386	1111.321	0.073	0.942
Allocationlithium:phasepost3	-3.749	4.259	1112.941	-0.880	0.379
Allocationlithium:phasepost4	-5.393	4.546	1118.762	-1.186	0.236
movement_typeL5:phasepre	1.522	3.704	1095.622	0.411	0.681
movement_typeL5:phasepost2	-4.326	4.204	1097.141	-1.029	0.304
movement_typeL5:phasepost3	-4.316	3.939	1097.240	-1.096	0.273
movement_typeL5:phasepost4	-6.258	4.297	1090.266	-1.456	0.146
Allocationlithium:movement_typeL5:phasepre	-7.817	5.474	1097.558	-1.428	0.154
Allocationlithium:movement_typeL5:phasepost2	0.374	6.075	1099.011	0.062	0.951
Allocationlithium:movement_typeL5:phasepost3	4.980	5.877	1098.714	0.847	0.397
Allocationlithium:movement_typeL5:phasepost4	7.112	6.221	1093.065	1.143	0.253

df_hlm_lmer_movement_mean$phase <- relevel(df_hlm_lmer_movement_mean$phase, ref = "post2")
kable(round(summary(lmer(model_formula_mean,data = df_hlm_lmer_movement_mean))$coefficients,3))

## boundary (singular) fit: see help('isSingular')

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	50.448	6.661	197.659	7.574	0.000
Allocationlithium	-8.934	6.097	51.213	-1.465	0.149
movement_typeL5	-59.777	5.041	67.253	-11.859	0.000
phasepost1	-3.444	3.024	1105.717	-1.139	0.255
phasepre	-5.807	2.867	1105.181	-2.026	0.043
phasepost3	-0.680	2.887	1097.231	-0.235	0.814
phasepost4	1.182	3.218	1111.685	0.367	0.713
posLog	4.494	1.458	757.526	3.082	0.002
pos_vmu	-0.613	0.454	804.907	-1.352	0.177
Age	-0.017	0.054	237.166	-0.312	0.756
GenderM	1.564	1.344	266.386	1.164	0.246
Allocationlithium:movement_typeL5	8.557	6.804	69.431	1.258	0.213
Allocationlithium:phasepost1	-0.320	4.386	1111.321	-0.073	0.942
Allocationlithium:phasepre	8.055	3.973	1106.937	2.028	0.043
Allocationlithium:phasepost3	-4.069	4.083	1099.931	-0.997	0.319
Allocationlithium:phasepost4	-5.713	4.480	1115.125	-1.275	0.203
movement_typeL5:phasepost1	4.326	4.204	1097.141	1.029	0.304
movement_typeL5:phasepre	5.849	3.970	1095.939	1.473	0.141
movement_typeL5:phasepost3	0.010	4.061	1095.973	0.002	0.998
movement_typeL5:phasepost4	-1.932	4.451	1093.847	-0.434	0.664
Allocationlithium:movement_typeL5:phasepost1	-0.374	6.075	1099.011	-0.062	0.951
Allocationlithium:movement_typeL5:phasepre	-8.191	5.525	1097.678	-1.482	0.139
Allocationlithium:movement_typeL5:phasepost3	4.606	5.726	1096.015	0.804	0.421
Allocationlithium:movement_typeL5:phasepost4	6.738	6.179	1094.876	1.091	0.276

df_hlm_lmer_movement_mean$phase <- relevel(df_hlm_lmer_movement_mean$phase, ref = "post3")
kable(round(summary(lmer(model_formula_mean,data = df_hlm_lmer_movement_mean))$coefficients,3))

## boundary (singular) fit: see help('isSingular')

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	49.769	6.579	189.266	7.565	0.000
Allocationlithium	-13.003	5.980	47.452	-2.174	0.035
movement_typeL5	-59.768	4.784	54.935	-12.493	0.000
phasepost2	0.680	2.887	1097.231	0.235	0.814
phasepost1	-2.765	2.851	1108.841	-0.970	0.332
phasepre	-5.127	2.652	1108.042	-1.933	0.053
phasepost4	1.862	2.978	1110.217	0.625	0.532
posLog	4.494	1.458	757.526	3.082	0.002
pos_vmu	-0.613	0.454	804.907	-1.352	0.177
Age	-0.017	0.054	237.166	-0.312	0.756
GenderM	1.564	1.344	266.386	1.164	0.246
Allocationlithium:movement_typeL5	13.163	6.598	61.663	1.995	0.050
Allocationlithium:phasepost2	4.069	4.083	1099.931	0.997	0.319
Allocationlithium:phasepost1	3.749	4.259	1112.941	0.880	0.379
Allocationlithium:phasepre	12.124	3.849	1111.723	3.150	0.002
Allocationlithium:phasepost4	-1.644	4.239	1111.287	-0.388	0.698
movement_typeL5:phasepost2	-0.010	4.061	1095.973	-0.002	0.998
movement_typeL5:phasepost1	4.316	3.939	1097.240	1.096	0.273
movement_typeL5:phasepre	5.839	3.670	1096.865	1.591	0.112
movement_typeL5:phasepost4	-1.942	4.139	1097.203	-0.469	0.639
Allocationlithium:movement_typeL5:phasepost2	-4.606	5.726	1096.015	-0.804	0.421
Allocationlithium:movement_typeL5:phasepost1	-4.980	5.877	1098.714	-0.847	0.397
Allocationlithium:movement_typeL5:phasepre	-12.797	5.315	1097.084	-2.408	0.016
Allocationlithium:movement_typeL5:phasepost4	2.132	5.903	1099.278	0.361	0.718

df_hlm_lmer_movement_mean$phase <- relevel(df_hlm_lmer_movement_mean$phase, ref = "post4")
kable(round(summary(lmer(model_formula_mean,data = df_hlm_lmer_movement_mean))$coefficients,3))

## boundary (singular) fit: see help('isSingular')

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	51.630	6.681	201.171	7.728	0.000
Allocationlithium	-14.647	6.129	52.572	-2.390	0.020
movement_typeL5	-61.710	5.035	67.324	-12.257	0.000
phasepost3	-1.862	2.978	1110.217	-0.625	0.532
phasepost2	-1.182	3.218	1111.685	-0.367	0.713
phasepost1	-4.626	3.138	1117.351	-1.474	0.141
phasepre	-6.989	2.963	1115.918	-2.359	0.019
posLog	4.494	1.458	757.526	3.082	0.002
pos_vmu	-0.613	0.454	804.907	-1.352	0.177
Age	-0.017	0.054	237.166	-0.312	0.756
GenderM	1.564	1.344	266.386	1.164	0.246
Allocationlithium:movement_typeL5	15.295	6.869	72.649	2.227	0.029
Allocationlithium:phasepost3	1.644	4.239	1111.287	0.388	0.698
Allocationlithium:phasepost2	5.713	4.480	1115.125	1.275	0.203
Allocationlithium:phasepost1	5.393	4.546	1118.762	1.186	0.236
Allocationlithium:phasepre	13.768	4.187	1119.807	3.288	0.001
movement_typeL5:phasepost3	1.942	4.139	1097.203	0.469	0.639
movement_typeL5:phasepost2	1.932	4.451	1093.847	0.434	0.664
movement_typeL5:phasepost1	6.258	4.297	1090.266	1.456	0.146
movement_typeL5:phasepre	7.781	4.053	1089.207	1.920	0.055
Allocationlithium:movement_typeL5:phasepost3	-2.132	5.903	1099.278	-0.361	0.718
Allocationlithium:movement_typeL5:phasepost2	-6.738	6.179	1094.876	-1.091	0.276
Allocationlithium:movement_typeL5:phasepost1	-7.112	6.221	1093.065	-1.143	0.253
Allocationlithium:movement_typeL5:phasepre	-14.929	5.720	1091.567	-2.610	0.009

df_hlm_lmer_movement_mean$phase <- factor(df_hlm_lmer_movement_mean$phase,
                                          levels = c("pre","post1","post2","post3","post4"))

plot_model(lmer(model_formula_mean,data = df_hlm_lmer_movement_mean),
           type = "pred",terms = c("phase","Allocation"))

## boundary (singular) fit: see help('isSingular')

# Plot Model using another method
df_hlm_lmer_movement_mean$movement_type <- relevel(df_hlm_lmer_movement_mean$movement_type, ref = "M10")
ee1 <- as.data.frame(Effect(c("phase", "Allocation","movement_type"),mod_movement_mean_M10))

stat.test_M10 <- as.data.frame(round(summary(mod_movement_mean_M10)$coefficients,3))[13:16,]
stat.test_L5 <- as.data.frame(round(summary(mod_movement_mean_L5)$coefficients,3))[13:16,]
stat.test <- rbind(stat.test_L5,stat.test_M10)%>%rename(p_val = 'Pr(>|t|)')
stat.test[1:4,"movement_type"] <- "L5" 
stat.test[5:8,"movement_type"] <- "M10" 

stat.test$group1 <- c("post1","post2","post3","post4")
stat.test$group2 <- c("post1","post2","post3","post4")

stat.test <- left_join(stat.test,ee1[ee1$phase != "pre" & ee1$Allocation == "control",]%>%
  rename(group1 = phase))%>%
  mutate(y.position = if_else(movement_type == "L5",ee1$upper[11],if_else(movement_type == "M10", ee1$upper[1],0)))%>%
  mutate(p_sig = if_else(p_val < 0.001,"***", if_else(p_val < 0.01, "**", if_else(p_val < 0.05, "*", "ns")) ))

## Joining, by = c("movement_type", "group1")

p4a <- ggplot(ee1, aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  facet_grid(~movement_type)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank(),
        axis.text.x = element_text(angle=45))+
  labs(title = "",y="movement level")+ 
  scale_y_continuous(limits = c(-10,80),breaks = seq(-10, 80, by = 10))+
  stat_pvalue_manual(stat.test, label = "p_sig",vjust = -2)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))
p4a

ggsave("../figs/fig2_hlm_movement_val.png",width = 6 ,height=4,dpi = 300)

p4a_M10 <- ggplot(ee1%>%filter(movement_type == "M10"), aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank(),
        axis.text.x = element_text(angle=45))+
  labs(title = "",y="daytime movement level", x ="")+ 
  scale_y_continuous(limits = c(40,85),breaks = seq(40, 85, by = 10))+
  stat_pvalue_manual(stat.test%>%filter(movement_type == "M10"), label = "p_sig",vjust = -4)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))
p4a_M10

variability

model_formula <- val ~ Allocation * var_type * movement_type * phase + posLog + pos_vmu + Age + Gender + (var_type * movement_type| ID)

df_hlm_lmer_movement_variability <- df_hlm%>%
  select(ID,Allocation,day,phase,movement_L5_s,movement_L5_vmu, movement_M10_s,movement_M10_vmu,Age,Gender,posLog, negLog, pos_vmu)%>%
  pivot_longer(c(movement_L5_s,movement_L5_vmu, movement_M10_s,movement_M10_vmu),names_pattern =  "movement_(.*)_(.*)",names_to = c("movement_type","var_type"),values_to = "val")%>%
  mutate(var_type = if_else(var_type == "vmu","volatility","noise"))

df_hlm_lmer_movement_variability$var_type <- factor(df_hlm_lmer_movement_variability$var_type)
df_hlm_lmer_movement_variability$var_type <- relevel(df_hlm_lmer_movement_variability$var_type, ref = "noise")
df_hlm_lmer_movement_variability$movement_type <- factor(df_hlm_lmer_movement_variability$movement_type)
df_hlm_lmer_movement_variability$movement_type <- relevel(df_hlm_lmer_movement_variability$movement_type, ref = "M10")

mod_movement_var_s_M10 <- lmer(model_formula,data = df_hlm_lmer_movement_variability)
kable(round(summary(mod_movement_var_s_M10)$coefficients,3))

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	-1.658	0.167	79.075	-9.919	0.000
Allocationlithium	-0.001	0.130	46.126	-0.009	0.993
var_typevolatility	-2.680	0.179	39.262	-15.002	0.000
movement_typeL5	0.078	0.137	45.907	0.571	0.571
phasepost1	-0.130	0.065	3379.204	-1.988	0.047
phasepost2	-0.090	0.067	3391.826	-1.348	0.178
phasepost3	-0.008	0.067	3388.243	-0.120	0.904
phasepost4	-0.009	0.068	3387.068	-0.134	0.894
posLog	0.097	0.028	3058.970	3.416	0.001
pos_vmu	0.028	0.009	2824.919	2.962	0.003
Age	0.000	0.003	28.449	-0.152	0.881
GenderM	0.127	0.073	30.311	1.743	0.091
Allocationlithium:var_typevolatility	0.013	0.240	40.105	0.054	0.957
Allocationlithium:movement_typeL5	-0.107	0.185	47.548	-0.580	0.565
var_typevolatility:movement_typeL5	-0.111	0.206	44.125	-0.537	0.594
Allocationlithium:phasepost1	0.136	0.090	3380.514	1.513	0.130
Allocationlithium:phasepost2	0.162	0.092	3383.483	1.755	0.079
Allocationlithium:phasepost3	-0.035	0.094	3387.121	-0.374	0.708
Allocationlithium:phasepost4	-0.082	0.096	3390.937	-0.854	0.393
var_typevolatility:phasepost1	-0.712	0.092	3374.134	-7.725	0.000
var_typevolatility:phasepost2	-0.922	0.094	3383.814	-9.773	0.000
var_typevolatility:phasepost3	-1.108	0.094	3381.921	-11.776	0.000
var_typevolatility:phasepost4	-1.242	0.096	3381.777	-12.929	0.000
movement_typeL5:phasepost1	-0.043	0.092	3375.088	-0.468	0.640
movement_typeL5:phasepost2	-0.055	0.094	3387.046	-0.582	0.561
movement_typeL5:phasepost3	-0.240	0.094	3384.797	-2.550	0.011
movement_typeL5:phasepost4	-0.264	0.096	3384.655	-2.746	0.006
Allocationlithium:var_typevolatility:movement_typeL5	-0.014	0.278	45.521	-0.049	0.961
Allocationlithium:var_typevolatility:phasepost1	0.314	0.127	3377.654	2.480	0.013
Allocationlithium:var_typevolatility:phasepost2	0.164	0.130	3379.757	1.259	0.208
Allocationlithium:var_typevolatility:phasepost3	0.264	0.133	3382.035	1.993	0.046
Allocationlithium:var_typevolatility:phasepost4	0.374	0.135	3385.154	2.777	0.006
Allocationlithium:movement_typeL5:phasepost1	-0.065	0.127	3379.591	-0.510	0.610
Allocationlithium:movement_typeL5:phasepost2	-0.076	0.130	3382.147	-0.579	0.562
Allocationlithium:movement_typeL5:phasepost3	0.179	0.133	3385.068	1.348	0.178
Allocationlithium:movement_typeL5:phasepost4	0.170	0.135	3388.911	1.263	0.207
var_typevolatility:movement_typeL5:phasepost1	0.184	0.130	3375.350	1.412	0.158
var_typevolatility:movement_typeL5:phasepost2	0.107	0.133	3387.487	0.806	0.421
var_typevolatility:movement_typeL5:phasepost3	0.353	0.133	3385.124	2.651	0.008
var_typevolatility:movement_typeL5:phasepost4	0.468	0.136	3384.987	3.447	0.001
Allocationlithium:var_typevolatility:movement_typeL5:phasepost1	-0.221	0.179	3379.929	-1.235	0.217
Allocationlithium:var_typevolatility:movement_typeL5:phasepost2	0.170	0.184	3382.511	0.923	0.356
Allocationlithium:var_typevolatility:movement_typeL5:phasepost3	-0.078	0.188	3385.482	-0.416	0.678
Allocationlithium:var_typevolatility:movement_typeL5:phasepost4	0.039	0.191	3389.365	0.203	0.839

df_hlm_lmer_movement_variability$var_type <- relevel(df_hlm_lmer_movement_variability$var_type, ref = "volatility")
df_hlm_lmer_movement_variability$movement_type <- relevel(df_hlm_lmer_movement_variability$movement_type, ref = "M10")

#same model, different reference level
mod_movement_var_vmu_M10 <- lmer(model_formula,data = df_hlm_lmer_movement_variability)
kable(round(summary(mod_movement_var_vmu_M10)$coefficients,3))

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	-4.338	0.216	73.002	-20.109	0.000
Allocationlithium	0.012	0.224	36.448	0.052	0.958
var_typenoise	2.680	0.179	39.265	15.002	0.000
movement_typeL5	-0.033	0.171	40.109	-0.191	0.850
phasepost1	-0.842	0.065	3375.844	-12.866	0.000
phasepost2	-1.013	0.067	3381.415	-15.079	0.000
phasepost3	-1.116	0.067	3379.165	-16.691	0.000
phasepost4	-1.251	0.068	3378.095	-18.353	0.000
posLog	0.097	0.028	3058.968	3.416	0.001
pos_vmu	0.028	0.009	2824.919	2.963	0.003
Age	0.000	0.003	28.447	-0.152	0.881
GenderM	0.127	0.073	30.309	1.743	0.091
Allocationlithium:var_typenoise	-0.013	0.240	40.109	-0.054	0.957
Allocationlithium:movement_typeL5	-0.121	0.230	41.052	-0.526	0.602
var_typenoise:movement_typeL5	0.111	0.206	44.121	0.537	0.594
Allocationlithium:phasepost1	0.450	0.090	3374.121	5.019	0.000
Allocationlithium:phasepost2	0.326	0.092	3375.347	3.533	0.000
Allocationlithium:phasepost3	0.229	0.094	3377.622	2.436	0.015
Allocationlithium:phasepost4	0.293	0.096	3379.557	3.059	0.002
var_typenoise:phasepost1	0.712	0.092	3374.131	7.725	0.000
var_typenoise:phasepost2	0.922	0.094	3383.811	9.773	0.000
var_typenoise:phasepost3	1.108	0.094	3381.918	11.776	0.000
var_typenoise:phasepost4	1.242	0.096	3381.774	12.929	0.000
movement_typeL5:phasepost1	0.141	0.092	3374.016	1.530	0.126
movement_typeL5:phasepost2	0.053	0.094	3383.395	0.557	0.578
movement_typeL5:phasepost3	0.113	0.094	3381.569	1.199	0.231
movement_typeL5:phasepost4	0.204	0.096	3381.452	2.127	0.034
Allocationlithium:var_typenoise:movement_typeL5	0.014	0.278	45.516	0.049	0.961
Allocationlithium:var_typenoise:phasepost1	-0.314	0.127	3377.651	-2.480	0.013
Allocationlithium:var_typenoise:phasepost2	-0.164	0.130	3379.754	-1.259	0.208
Allocationlithium:var_typenoise:phasepost3	-0.264	0.133	3382.032	-1.993	0.046
Allocationlithium:var_typenoise:phasepost4	-0.374	0.135	3385.151	-2.777	0.006
Allocationlithium:movement_typeL5:phasepost1	-0.286	0.127	3377.448	-2.257	0.024
Allocationlithium:movement_typeL5:phasepost2	0.095	0.130	3379.453	0.725	0.468
Allocationlithium:movement_typeL5:phasepost3	0.101	0.133	3381.764	0.760	0.447
Allocationlithium:movement_typeL5:phasepost4	0.209	0.135	3384.917	1.550	0.121
var_typenoise:movement_typeL5:phasepost1	-0.184	0.130	3375.347	-1.412	0.158
var_typenoise:movement_typeL5:phasepost2	-0.107	0.133	3387.484	-0.806	0.421
var_typenoise:movement_typeL5:phasepost3	-0.353	0.133	3385.121	-2.651	0.008
var_typenoise:movement_typeL5:phasepost4	-0.468	0.136	3384.984	-3.447	0.001
Allocationlithium:var_typenoise:movement_typeL5:phasepost1	0.221	0.179	3379.926	1.235	0.217
Allocationlithium:var_typenoise:movement_typeL5:phasepost2	-0.170	0.184	3382.509	-0.923	0.356
Allocationlithium:var_typenoise:movement_typeL5:phasepost3	0.078	0.188	3385.480	0.416	0.678
Allocationlithium:var_typenoise:movement_typeL5:phasepost4	-0.039	0.191	3389.363	-0.203	0.839

df_hlm_lmer_movement_variability$Allocation <- factor(df_hlm_lmer_movement_variability$Allocation, levels = c("lithium","control"))
kable(round(summary(lmer(model_formula,data = df_hlm_lmer_movement_variability))$coefficients,3))

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	-4.326	0.195	80.657	-22.234	0.000
Allocationcontrol	-0.012	0.224	36.461	-0.052	0.958
var_typenoise	2.667	0.161	41.203	16.600	0.000
movement_typeL5	-0.153	0.154	42.279	-0.998	0.324
phasepost1	-0.392	0.062	3378.474	-6.361	0.000
phasepost2	-0.687	0.064	3376.060	-10.730	0.000
phasepost3	-0.887	0.066	3376.527	-13.388	0.000
phasepost4	-0.959	0.067	3380.460	-14.305	0.000
posLog	0.097	0.028	3058.989	3.416	0.001
pos_vmu	0.028	0.009	2824.958	2.963	0.003
Age	0.000	0.003	28.447	-0.152	0.880
GenderM	0.127	0.073	30.310	1.743	0.091
Allocationcontrol:var_typenoise	0.013	0.240	40.121	0.054	0.957
Allocationcontrol:movement_typeL5	0.121	0.230	41.069	0.526	0.602
var_typenoise:movement_typeL5	0.124	0.186	47.331	0.667	0.508
Allocationcontrol:phasepost1	-0.450	0.090	3374.114	-5.018	0.000
Allocationcontrol:phasepost2	-0.326	0.092	3375.340	-3.533	0.000
Allocationcontrol:phasepost3	-0.229	0.094	3377.615	-2.436	0.015
Allocationcontrol:phasepost4	-0.293	0.096	3379.551	-3.059	0.002
var_typenoise:phasepost1	0.398	0.087	3381.330	4.577	0.000
var_typenoise:phasepost2	0.758	0.090	3374.860	8.416	0.000
var_typenoise:phasepost3	0.844	0.094	3382.309	9.014	0.000
var_typenoise:phasepost4	0.868	0.095	3388.568	9.178	0.000
movement_typeL5:phasepost1	-0.145	0.087	3381.050	-1.667	0.096
movement_typeL5:phasepost2	0.147	0.090	3374.709	1.634	0.102
movement_typeL5:phasepost3	0.214	0.094	3381.972	2.283	0.022
movement_typeL5:phasepost4	0.413	0.095	3388.198	4.370	0.000
Allocationcontrol:var_typenoise:movement_typeL5	-0.014	0.278	45.534	-0.049	0.961
Allocationcontrol:var_typenoise:phasepost1	0.314	0.127	3377.644	2.480	0.013
Allocationcontrol:var_typenoise:phasepost2	0.164	0.130	3379.748	1.259	0.208
Allocationcontrol:var_typenoise:phasepost3	0.264	0.133	3382.026	1.993	0.046
Allocationcontrol:var_typenoise:phasepost4	0.374	0.135	3385.145	2.777	0.006
Allocationcontrol:movement_typeL5:phasepost1	0.286	0.127	3377.442	2.257	0.024
Allocationcontrol:movement_typeL5:phasepost2	-0.095	0.130	3379.447	-0.725	0.468
Allocationcontrol:movement_typeL5:phasepost3	-0.101	0.133	3381.757	-0.760	0.447
Allocationcontrol:movement_typeL5:phasepost4	-0.209	0.135	3384.911	-1.550	0.121
var_typenoise:movement_typeL5:phasepost1	0.037	0.123	3384.588	0.303	0.762
var_typenoise:movement_typeL5:phasepost2	-0.278	0.127	3376.255	-2.180	0.029
var_typenoise:movement_typeL5:phasepost3	-0.275	0.132	3385.856	-2.077	0.038
var_typenoise:movement_typeL5:phasepost4	-0.507	0.134	3393.269	-3.792	0.000
Allocationcontrol:var_typenoise:movement_typeL5:phasepost1	-0.221	0.179	3379.920	-1.235	0.217
Allocationcontrol:var_typenoise:movement_typeL5:phasepost2	0.170	0.184	3382.502	0.923	0.356
Allocationcontrol:var_typenoise:movement_typeL5:phasepost3	-0.078	0.188	3385.474	-0.416	0.678
Allocationcontrol:var_typenoise:movement_typeL5:phasepost4	0.039	0.191	3389.357	0.203	0.839

df_hlm_lmer_movement_variability$Allocation <- relevel(df_hlm_lmer_movement_variability$Allocation, ref = "control")

#same model, different reference level
df_hlm_lmer_movement_variability$var_type <- relevel(df_hlm_lmer_movement_variability$var_type, ref = "noise")
df_hlm_lmer_movement_variability$movement_type <- relevel(df_hlm_lmer_movement_variability$movement_type, ref = "L5")

mod_movement_var_s_L5 <- lmer(model_formula,data = df_hlm_lmer_movement_variability)
kable(round(summary(mod_movement_var_s_L5)$coefficients,3))

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	-1.580	0.163	75.855	-9.703	0.000
Allocationlithium	-0.108	0.120	49.484	-0.902	0.371
var_typevolatility	-2.791	0.122	49.626	-22.859	0.000
movement_typeM10	-0.078	0.137	45.901	-0.571	0.571
phasepost1	-0.173	0.065	3381.659	-2.647	0.008
phasepost2	-0.145	0.067	3397.093	-2.168	0.030
phasepost3	-0.248	0.067	3393.075	-3.713	0.000
phasepost4	-0.273	0.068	3391.963	-4.009	0.000
posLog	0.097	0.028	3058.974	3.416	0.001
pos_vmu	0.028	0.009	2824.919	2.963	0.003
Age	0.000	0.003	28.450	-0.152	0.880
GenderM	0.127	0.073	30.313	1.743	0.091
Allocationlithium:var_typevolatility	-0.001	0.165	51.772	-0.004	0.997
Allocationlithium:movement_typeM10	0.107	0.185	47.542	0.580	0.565
var_typevolatility:movement_typeM10	0.111	0.206	44.117	0.537	0.594
Allocationlithium:phasepost1	0.071	0.090	3384.558	0.792	0.429
Allocationlithium:phasepost2	0.086	0.092	3388.230	0.937	0.349
Allocationlithium:phasepost3	0.144	0.094	3392.444	1.528	0.127
Allocationlithium:phasepost4	0.089	0.095	3396.827	0.928	0.353
var_typevolatility:phasepost1	-0.528	0.092	3378.822	-5.730	0.000
var_typevolatility:phasepost2	-0.815	0.094	3395.366	-8.647	0.000
var_typevolatility:phasepost3	-0.755	0.094	3392.150	-8.038	0.000
var_typevolatility:phasepost4	-0.774	0.096	3391.974	-8.067	0.000
movement_typeM10:phasepost1	0.043	0.092	3375.091	0.468	0.640
movement_typeM10:phasepost2	0.055	0.094	3387.049	0.582	0.561
movement_typeM10:phasepost3	0.240	0.094	3384.800	2.550	0.011
movement_typeM10:phasepost4	0.264	0.096	3384.657	2.746	0.006
Allocationlithium:var_typevolatility:movement_typeM10	0.014	0.278	45.512	0.049	0.961
Allocationlithium:var_typevolatility:phasepost1	0.093	0.127	3385.508	0.734	0.463
Allocationlithium:var_typevolatility:phasepost2	0.335	0.130	3389.064	2.566	0.010
Allocationlithium:var_typevolatility:phasepost3	0.187	0.133	3392.905	1.407	0.160
Allocationlithium:var_typevolatility:phasepost4	0.413	0.135	3397.488	3.069	0.002
Allocationlithium:movement_typeM10:phasepost1	0.065	0.127	3379.594	0.510	0.610
Allocationlithium:movement_typeM10:phasepost2	0.076	0.130	3382.150	0.579	0.562
Allocationlithium:movement_typeM10:phasepost3	-0.179	0.133	3385.071	-1.348	0.178
Allocationlithium:movement_typeM10:phasepost4	-0.170	0.135	3388.913	-1.263	0.207
var_typevolatility:movement_typeM10:phasepost1	-0.184	0.130	3375.353	-1.412	0.158
var_typevolatility:movement_typeM10:phasepost2	-0.107	0.133	3387.489	-0.806	0.421
var_typevolatility:movement_typeM10:phasepost3	-0.353	0.133	3385.126	-2.651	0.008
var_typevolatility:movement_typeM10:phasepost4	-0.468	0.136	3384.989	-3.447	0.001
Allocationlithium:var_typevolatility:movement_typeM10:phasepost1	0.221	0.179	3379.932	1.235	0.217
Allocationlithium:var_typevolatility:movement_typeM10:phasepost2	-0.170	0.184	3382.514	-0.923	0.356
Allocationlithium:var_typevolatility:movement_typeM10:phasepost3	0.078	0.188	3385.485	0.416	0.678
Allocationlithium:var_typevolatility:movement_typeM10:phasepost4	-0.039	0.191	3389.367	-0.203	0.839

df_hlm_lmer_movement_variability$var_type <- relevel(df_hlm_lmer_movement_variability$var_type, ref = "volatility")
df_hlm_lmer_movement_variability$movement_type <- relevel(df_hlm_lmer_movement_variability$movement_type, ref = "L5")

mod_movement_var_vmu_L5 <- lmer(model_formula,data = df_hlm_lmer_movement_variability)
kable(round(summary(mod_movement_var_vmu_L5)$coefficients,3))

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	-4.370	0.188	82.911	-23.216	0.000
Allocationlithium	-0.109	0.175	39.862	-0.624	0.537
var_typenoise	2.791	0.122	49.639	22.861	0.000
movement_typeM10	0.033	0.171	40.112	0.191	0.850
phasepost1	-0.701	0.065	3377.827	-10.713	0.000
phasepost2	-0.960	0.067	3388.169	-14.305	0.000
phasepost3	-1.003	0.067	3384.895	-15.012	0.000
phasepost4	-1.047	0.068	3383.872	-15.365	0.000
posLog	0.097	0.028	3058.972	3.416	0.001
pos_vmu	0.028	0.009	2824.920	2.963	0.003
Age	0.000	0.003	28.449	-0.152	0.880
GenderM	0.127	0.073	30.312	1.743	0.091
Allocationlithium:var_typenoise	0.001	0.165	51.786	0.004	0.997
Allocationlithium:movement_typeM10	0.121	0.230	41.055	0.526	0.602
var_typenoise:movement_typeM10	-0.111	0.206	44.124	-0.537	0.594
Allocationlithium:phasepost1	0.164	0.090	3378.038	1.829	0.068
Allocationlithium:phasepost2	0.421	0.092	3380.176	4.559	0.000
Allocationlithium:phasepost3	0.330	0.094	3383.670	3.509	0.000
Allocationlithium:phasepost4	0.502	0.096	3387.088	5.246	0.000
var_typenoise:phasepost1	0.528	0.092	3378.819	5.730	0.000
var_typenoise:phasepost2	0.815	0.094	3395.363	8.647	0.000
var_typenoise:phasepost3	0.755	0.094	3392.147	8.038	0.000
var_typenoise:phasepost4	0.774	0.096	3391.971	8.067	0.000
movement_typeM10:phasepost1	-0.141	0.092	3374.017	-1.529	0.126
movement_typeM10:phasepost2	-0.053	0.094	3383.396	-0.557	0.578
movement_typeM10:phasepost3	-0.113	0.094	3381.571	-1.199	0.231
movement_typeM10:phasepost4	-0.204	0.096	3381.453	-2.127	0.034
Allocationlithium:var_typenoise:movement_typeM10	-0.014	0.278	45.520	-0.049	0.961
Allocationlithium:var_typenoise:phasepost1	-0.093	0.127	3385.505	-0.734	0.463
Allocationlithium:var_typenoise:phasepost2	-0.335	0.130	3389.061	-2.566	0.010
Allocationlithium:var_typenoise:phasepost3	-0.187	0.133	3392.902	-1.407	0.160
Allocationlithium:var_typenoise:phasepost4	-0.413	0.135	3397.485	-3.069	0.002
Allocationlithium:movement_typeM10:phasepost1	0.286	0.127	3377.449	2.257	0.024
Allocationlithium:movement_typeM10:phasepost2	-0.095	0.130	3379.455	-0.725	0.468
Allocationlithium:movement_typeM10:phasepost3	-0.101	0.133	3381.765	-0.760	0.447
Allocationlithium:movement_typeM10:phasepost4	-0.209	0.135	3384.919	-1.550	0.121
var_typenoise:movement_typeM10:phasepost1	0.184	0.130	3375.349	1.412	0.158
var_typenoise:movement_typeM10:phasepost2	0.107	0.133	3387.486	0.806	0.421
var_typenoise:movement_typeM10:phasepost3	0.353	0.133	3385.123	2.651	0.008
var_typenoise:movement_typeM10:phasepost4	0.468	0.136	3384.986	3.447	0.001
Allocationlithium:var_typenoise:movement_typeM10:phasepost1	-0.221	0.179	3379.928	-1.235	0.217
Allocationlithium:var_typenoise:movement_typeM10:phasepost2	0.170	0.184	3382.511	0.923	0.356
Allocationlithium:var_typenoise:movement_typeM10:phasepost3	-0.078	0.188	3385.482	-0.416	0.678
Allocationlithium:var_typenoise:movement_typeM10:phasepost4	0.039	0.191	3389.365	0.203	0.839

plot_model(mod_movement_var_s_M10, type = "pred", terms = c("phase", "Allocation","movement_type","var_type"))

Conclusion: lithium increased volatility of M10 and L5 by week 4. Conclusion: Lithium increases daytime volatility and decreases nighttime volatility of time

relevel to different weeks to show group contrast

df_hlm_lmer_movement_variability$movement_type <- relevel(df_hlm_lmer_movement_variability$movement_type, ref = "M10")
df_hlm_lmer_movement_variability$phase <- relevel(df_hlm_lmer_movement_variability$phase, ref = "post1")
kable(round(summary(lmer(model_formula,data = df_hlm_lmer_movement_variability))$coefficients,3))

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	-5.180	0.216	73.680	-23.959	0.000
Allocationlithium	0.462	0.225	36.779	2.052	0.047
var_typenoise	3.392	0.180	40.332	18.863	0.000
movement_typeL5	0.108	0.172	41.306	0.631	0.532
phasepre	0.842	0.065	3375.844	12.866	0.000
phasepost2	-0.171	0.068	3376.764	-2.501	0.012
phasepost3	-0.274	0.068	3376.823	-4.014	0.000
phasepost4	-0.409	0.070	3377.550	-5.858	0.000
posLog	0.097	0.028	3058.986	3.416	0.001
pos_vmu	0.028	0.009	2824.942	2.962	0.003
Age	0.000	0.003	28.451	-0.152	0.881
GenderM	0.127	0.073	30.313	1.743	0.091
Allocationlithium:var_typenoise	-0.327	0.241	40.737	-1.357	0.182
Allocationlithium:movement_typeL5	-0.407	0.231	41.759	-1.764	0.085
var_typenoise:movement_typeL5	-0.073	0.208	45.936	-0.353	0.726
Allocationlithium:phasepre	-0.450	0.090	3374.121	-5.019	0.000
Allocationlithium:phasepost2	-0.123	0.094	3375.672	-1.314	0.189
Allocationlithium:phasepost3	-0.220	0.096	3376.613	-2.307	0.021
Allocationlithium:phasepost4	-0.157	0.097	3376.411	-1.622	0.105
var_typenoise:phasepre	-0.712	0.092	3374.131	-7.725	0.000
var_typenoise:phasepost2	0.210	0.096	3381.885	2.183	0.029
var_typenoise:phasepost3	0.396	0.096	3381.518	4.111	0.000
var_typenoise:phasepost4	0.530	0.099	3382.890	5.380	0.000
movement_typeL5:phasepre	-0.141	0.092	3374.016	-1.529	0.126
movement_typeL5:phasepost2	-0.088	0.096	3381.553	-0.917	0.359
movement_typeL5:phasepost3	-0.028	0.096	3381.224	-0.292	0.770
movement_typeL5:phasepost4	0.063	0.099	3382.585	0.643	0.520
Allocationlithium:var_typenoise:movement_typeL5	0.235	0.279	46.623	0.841	0.405
Allocationlithium:var_typenoise:phasepre	0.314	0.127	3377.651	2.480	0.013
Allocationlithium:var_typenoise:phasepost2	0.150	0.133	3380.229	1.129	0.259
Allocationlithium:var_typenoise:phasepost3	0.050	0.135	3380.805	0.369	0.712
Allocationlithium:var_typenoise:phasepost4	-0.060	0.136	3380.010	-0.440	0.660
Allocationlithium:movement_typeL5:phasepre	0.286	0.127	3377.448	2.257	0.024
Allocationlithium:movement_typeL5:phasepost2	0.381	0.133	3379.933	2.866	0.004
Allocationlithium:movement_typeL5:phasepost3	0.387	0.135	3380.557	2.869	0.004
Allocationlithium:movement_typeL5:phasepost4	0.495	0.136	3379.924	3.626	0.000
var_typenoise:movement_typeL5:phasepre	0.184	0.130	3375.347	1.412	0.158
var_typenoise:movement_typeL5:phasepost2	0.077	0.136	3385.013	0.562	0.574
var_typenoise:movement_typeL5:phasepost3	-0.169	0.136	3384.586	-1.238	0.216
var_typenoise:movement_typeL5:phasepost4	-0.284	0.139	3386.394	-2.039	0.042
Allocationlithium:var_typenoise:movement_typeL5:phasepre	-0.221	0.179	3379.927	-1.235	0.217
Allocationlithium:var_typenoise:movement_typeL5:phasepost2	-0.391	0.188	3383.024	-2.086	0.037
Allocationlithium:var_typenoise:movement_typeL5:phasepost3	-0.143	0.191	3383.830	-0.752	0.452
Allocationlithium:var_typenoise:movement_typeL5:phasepost4	-0.260	0.193	3383.081	-1.348	0.178

df_hlm_lmer_movement_variability$phase <- relevel(df_hlm_lmer_movement_variability$phase, ref = "post2")
kable(round(summary(lmer(model_formula,data = df_hlm_lmer_movement_variability))$coefficients,3))

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	-5.350	0.218	75.992	-24.553	0.000
Allocationlithium	0.338	0.226	37.652	1.494	0.144
var_typenoise	3.602	0.181	41.396	19.897	0.000
movement_typeL5	0.020	0.173	42.502	0.116	0.908
phasepost1	0.171	0.068	3376.764	2.501	0.012
phasepre	1.013	0.067	3381.416	15.078	0.000
phasepost3	-0.103	0.069	3371.569	-1.504	0.133
phasepost4	-0.238	0.070	3373.140	-3.389	0.001
posLog	0.097	0.028	3058.975	3.416	0.001
pos_vmu	0.028	0.009	2824.933	2.963	0.003
Age	0.000	0.003	28.446	-0.152	0.880
GenderM	0.127	0.073	30.308	1.743	0.091
Allocationlithium:var_typenoise	-0.177	0.244	42.429	-0.728	0.471
Allocationlithium:movement_typeL5	-0.026	0.233	43.661	-0.112	0.912
var_typenoise:movement_typeL5	0.003	0.210	47.683	0.015	0.988
Allocationlithium:phasepost1	0.123	0.094	3375.673	1.314	0.189
Allocationlithium:phasepre	-0.326	0.092	3375.348	-3.533	0.000
Allocationlithium:phasepost3	-0.097	0.098	3373.944	-0.994	0.320
Allocationlithium:phasepost4	-0.034	0.099	3375.862	-0.338	0.735
var_typenoise:phasepost1	-0.210	0.096	3381.886	-2.183	0.029
var_typenoise:phasepre	-0.922	0.094	3383.812	-9.773	0.000
var_typenoise:phasepost3	0.186	0.097	3372.000	1.912	0.056
var_typenoise:phasepost4	0.320	0.099	3374.944	3.216	0.001
movement_typeL5:phasepost1	0.088	0.096	3381.554	0.917	0.359
movement_typeL5:phasepre	-0.053	0.094	3383.396	-0.557	0.578
movement_typeL5:phasepost3	0.060	0.097	3371.959	0.621	0.535
movement_typeL5:phasepost4	0.152	0.099	3374.797	1.527	0.127
Allocationlithium:var_typenoise:movement_typeL5	-0.157	0.284	49.457	-0.552	0.583
Allocationlithium:var_typenoise:phasepost1	-0.150	0.133	3380.230	-1.129	0.259
Allocationlithium:var_typenoise:phasepre	0.164	0.130	3379.755	1.259	0.208
Allocationlithium:var_typenoise:phasepost3	-0.100	0.138	3376.294	-0.727	0.467
Allocationlithium:var_typenoise:phasepost4	-0.210	0.140	3379.967	-1.502	0.133
Allocationlithium:movement_typeL5:phasepost1	-0.381	0.133	3379.934	-2.866	0.004
Allocationlithium:movement_typeL5:phasepre	-0.095	0.130	3379.454	-0.725	0.468
Allocationlithium:movement_typeL5:phasepost3	0.006	0.138	3376.122	0.045	0.964
Allocationlithium:movement_typeL5:phasepost4	0.114	0.140	3379.764	0.817	0.414
var_typenoise:movement_typeL5:phasepost1	-0.077	0.136	3385.014	-0.562	0.574
var_typenoise:movement_typeL5:phasepre	0.107	0.133	3387.486	0.806	0.421
var_typenoise:movement_typeL5:phasepost3	-0.245	0.137	3372.430	-1.786	0.074
var_typenoise:movement_typeL5:phasepost4	-0.361	0.141	3376.433	-2.566	0.010
Allocationlithium:var_typenoise:movement_typeL5:phasepost1	0.391	0.188	3383.025	2.086	0.037
Allocationlithium:var_typenoise:movement_typeL5:phasepre	0.170	0.184	3382.510	0.923	0.356
Allocationlithium:var_typenoise:movement_typeL5:phasepost3	0.248	0.195	3378.150	1.274	0.203
Allocationlithium:var_typenoise:movement_typeL5:phasepost4	0.132	0.198	3382.951	0.665	0.506

df_hlm_lmer_movement_variability$phase <- relevel(df_hlm_lmer_movement_variability$phase, ref = "post3")
kable(round(summary(lmer(model_formula,data = df_hlm_lmer_movement_variability))$coefficients,3))

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00266921 (tol = 0.002, component 1)

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	-5.454	0.218	76.348	-25.001	0.000
Allocationlithium	0.241	0.227	38.148	1.062	0.295
var_typenoise	3.788	0.181	41.355	20.929	0.000
movement_typeL5	0.080	0.173	42.452	0.464	0.645
phasepost2	0.103	0.069	3371.568	1.504	0.133
phasepost1	0.274	0.068	3376.823	4.014	0.000
phasepre	1.116	0.067	3379.165	16.691	0.000
phasepost4	-0.135	0.070	3371.847	-1.926	0.054
posLog	0.097	0.028	3058.986	3.416	0.001
pos_vmu	0.028	0.009	2824.941	2.962	0.003
Age	0.000	0.003	28.450	-0.152	0.881
GenderM	0.127	0.073	30.312	1.743	0.091
Allocationlithium:var_typenoise	-0.277	0.245	43.274	-1.133	0.263
Allocationlithium:movement_typeL5	-0.020	0.234	44.607	-0.085	0.933
var_typenoise:movement_typeL5	-0.242	0.210	47.611	-1.152	0.255
Allocationlithium:phasepost2	0.097	0.098	3373.943	0.994	0.320
Allocationlithium:phasepost1	0.220	0.096	3376.613	2.307	0.021
Allocationlithium:phasepre	-0.229	0.094	3377.622	-2.436	0.015
Allocationlithium:phasepost4	0.063	0.100	3373.835	0.635	0.525
var_typenoise:phasepost2	-0.186	0.097	3372.000	-1.912	0.056
var_typenoise:phasepost1	-0.396	0.096	3381.518	-4.111	0.000
var_typenoise:phasepre	-1.108	0.094	3381.919	-11.776	0.000
var_typenoise:phasepost4	0.134	0.099	3372.954	1.352	0.177
movement_typeL5:phasepost2	-0.060	0.097	3371.958	-0.621	0.535
movement_typeL5:phasepost1	0.028	0.096	3381.224	0.292	0.770
movement_typeL5:phasepre	-0.113	0.094	3381.570	-1.199	0.231
movement_typeL5:phasepost4	0.091	0.099	3372.877	0.923	0.356
Allocationlithium:var_typenoise:movement_typeL5	0.092	0.286	50.888	0.321	0.750
Allocationlithium:var_typenoise:phasepost2	0.100	0.138	3376.293	0.727	0.467
Allocationlithium:var_typenoise:phasepost1	-0.050	0.135	3380.805	-0.369	0.712
Allocationlithium:var_typenoise:phasepre	0.264	0.133	3382.032	1.993	0.046
Allocationlithium:var_typenoise:phasepost4	-0.110	0.141	3376.995	-0.778	0.437
Allocationlithium:movement_typeL5:phasepost2	-0.006	0.138	3376.121	-0.045	0.964
Allocationlithium:movement_typeL5:phasepost1	-0.387	0.135	3380.558	-2.869	0.004
Allocationlithium:movement_typeL5:phasepre	-0.101	0.133	3381.764	-0.760	0.447
Allocationlithium:movement_typeL5:phasepost4	0.108	0.141	3376.777	0.765	0.444
var_typenoise:movement_typeL5:phasepost2	0.245	0.137	3372.430	1.786	0.074
var_typenoise:movement_typeL5:phasepost1	0.169	0.136	3384.587	1.238	0.216
var_typenoise:movement_typeL5:phasepre	0.353	0.133	3385.122	2.651	0.008
var_typenoise:movement_typeL5:phasepost4	-0.115	0.140	3373.736	-0.823	0.411
Allocationlithium:var_typenoise:movement_typeL5:phasepost2	-0.248	0.195	3378.149	-1.274	0.203
Allocationlithium:var_typenoise:movement_typeL5:phasepost1	0.143	0.191	3383.831	0.752	0.452
Allocationlithium:var_typenoise:movement_typeL5:phasepre	-0.078	0.188	3385.480	-0.416	0.678
Allocationlithium:var_typenoise:movement_typeL5:phasepost4	-0.117	0.200	3378.935	-0.584	0.559

df_hlm_lmer_movement_variability$phase <- relevel(df_hlm_lmer_movement_variability$phase, ref = "post4")
kable(round(summary(lmer(model_formula,data = df_hlm_lmer_movement_variability))$coefficients,3))

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00341248 (tol = 0.002, component 1)

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	-5.589	0.219	76.641	-25.568	0.000
Allocationlithium	0.304	0.228	38.323	1.338	0.189
var_typenoise	3.922	0.182	42.263	21.525	0.000
movement_typeL5	0.172	0.174	43.466	0.985	0.330
phasepost3	0.135	0.070	3371.897	1.927	0.054
phasepost2	0.238	0.070	3373.189	3.389	0.001
phasepost1	0.409	0.070	3377.596	5.858	0.000
phasepre	1.251	0.068	3378.141	18.353	0.000
posLog	0.097	0.028	3059.012	3.416	0.001
pos_vmu	0.028	0.009	2824.966	2.963	0.003
Age	0.000	0.003	28.435	-0.152	0.880
GenderM	0.127	0.073	30.296	1.743	0.091
Allocationlithium:var_typenoise	-0.387	0.246	43.683	-1.576	0.122
Allocationlithium:movement_typeL5	0.088	0.235	45.061	0.375	0.710
var_typenoise:movement_typeL5	-0.357	0.212	49.208	-1.685	0.098
Allocationlithium:phasepost3	-0.063	0.100	3373.883	-0.635	0.525
Allocationlithium:phasepost2	0.034	0.099	3375.909	0.338	0.735
Allocationlithium:phasepost1	0.157	0.097	3376.459	1.622	0.105
Allocationlithium:phasepre	-0.293	0.096	3379.602	-3.059	0.002
var_typenoise:phasepost3	-0.134	0.099	3373.004	-1.352	0.177
var_typenoise:phasepost2	-0.320	0.099	3374.993	-3.216	0.001
var_typenoise:phasepost1	-0.530	0.099	3382.936	-5.380	0.000
var_typenoise:phasepre	-1.242	0.096	3381.821	-12.929	0.000
movement_typeL5:phasepost3	-0.091	0.099	3372.926	-0.923	0.356
movement_typeL5:phasepost2	-0.152	0.099	3374.845	-1.527	0.127
movement_typeL5:phasepost1	-0.063	0.099	3382.630	-0.643	0.520
movement_typeL5:phasepre	-0.204	0.096	3381.497	-2.127	0.034
Allocationlithium:var_typenoise:movement_typeL5	-0.025	0.287	51.668	-0.087	0.931
Allocationlithium:var_typenoise:phasepost3	0.110	0.141	3377.044	0.778	0.437
Allocationlithium:var_typenoise:phasepost2	0.210	0.140	3380.014	1.502	0.133
Allocationlithium:var_typenoise:phasepost1	0.060	0.136	3380.058	0.440	0.660
Allocationlithium:var_typenoise:phasepre	0.374	0.135	3385.197	2.777	0.006
Allocationlithium:movement_typeL5:phasepost3	-0.108	0.141	3376.824	-0.765	0.444
Allocationlithium:movement_typeL5:phasepost2	-0.114	0.140	3379.809	-0.817	0.414
Allocationlithium:movement_typeL5:phasepost1	-0.495	0.136	3379.970	-3.626	0.000
Allocationlithium:movement_typeL5:phasepre	-0.209	0.135	3384.960	-1.550	0.121
var_typenoise:movement_typeL5:phasepost3	0.115	0.140	3373.786	0.823	0.411
var_typenoise:movement_typeL5:phasepost2	0.361	0.141	3376.481	2.566	0.010
var_typenoise:movement_typeL5:phasepost1	0.284	0.139	3386.439	2.039	0.042
var_typenoise:movement_typeL5:phasepre	0.468	0.136	3385.030	3.447	0.001
Allocationlithium:var_typenoise:movement_typeL5:phasepost3	0.117	0.200	3378.983	0.584	0.559
Allocationlithium:var_typenoise:movement_typeL5:phasepost2	-0.132	0.198	3382.997	-0.665	0.506
Allocationlithium:var_typenoise:movement_typeL5:phasepost1	0.260	0.193	3383.127	1.348	0.178
Allocationlithium:var_typenoise:movement_typeL5:phasepre	0.039	0.191	3389.407	0.203	0.839

df_hlm_lmer_movement_variability$phase <- factor(df_hlm_lmer_movement_variability$phase,
                                          levels = c("pre","post1","post2","post3","post4"))

plot_model(lmer(model_formula,data = df_hlm_lmer_movement_variability),
           type = "pred",terms = c("phase","Allocation"))

df_hlm_lmer_movement_variability$var_type <- relevel(df_hlm_lmer_movement_variability$var_type, ref = "volatility")
df_hlm_lmer_movement_variability$movement_type <- relevel(df_hlm_lmer_movement_variability$movement_type, ref = "M10")
ee2 <- as.data.frame(Effect(c("phase", "Allocation","movement_type","var_type"),mod_movement_var_vmu_M10))

stat.test_s_M10 <- as.data.frame(round(summary(mod_movement_var_s_M10)$coefficients,3))[16:19,]
stat.test_s_L5 <- as.data.frame(round(summary(mod_movement_var_s_L5)$coefficients,3))[16:19,]
stat.test_vmu_M10 <- as.data.frame(round(summary(mod_movement_var_vmu_M10)$coefficients,3))[16:19,]
stat.test_vmu_L5 <- as.data.frame(round(summary(mod_movement_var_vmu_L5)$coefficients,3))[16:19,]

stat.test <- Reduce(rbind,list(stat.test_s_M10,stat.test_s_L5,stat.test_vmu_M10,stat.test_vmu_L5))%>%rename(p_val = 'Pr(>|t|)')
stat.test[c(1:4,9:12),"movement_type"] <- "M10" 
stat.test[c(5:8,13:16),"movement_type"] <- "L5" 

stat.test[1:8,"var_type"] <- "noise" 
stat.test[9:16,"var_type"] <- "volatility" 

stat.test$group1 <- c("post1","post2","post3","post4")
stat.test$group2 <- c("post1","post2","post3","post4")

stat.test <- left_join(stat.test,ee2[ee2$phase != "pre" & ee2$Allocation == "lithium",]%>%
  rename(group1 = phase), by = c("movement_type","var_type","group1"))%>%
  mutate(y.position = if_else(var_type == "volatility",ee2$upper[1],if_else(var_type == "noise", ee2$upper[21],0)))%>%
  mutate(p_sig = if_else(p_val < 0.001,"***", if_else(p_val < 0.01, "**", if_else(p_val < 0.05, "*", "ns")) ))

#plot noise first
p11a <- ggplot(ee2%>%filter(var_type == "noise"), aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  facet_wrap(~movement_type)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank())+
  labs(title = "",y="movement noise")+
  scale_y_continuous(limits = c(-2,-1),breaks = seq(-2,-1, by = 0.2))+
  stat_pvalue_manual(stat.test%>%filter(var_type == "noise"),label = "p_sig",vjust = -3)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))


#plot volatility first
p11b <- ggplot(ee2%>%filter(var_type == "volatility"), aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  facet_wrap(~movement_type)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank())+
  labs(title = "",y="movement volatility")+
  scale_y_continuous(limits = c(-5.8,-3.8),breaks = seq(-5.8,-3.8, by = 0.4))+
  stat_pvalue_manual(stat.test%>%filter(var_type == "volatility"),label = "p_sig",vjust = 1)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))

p11a/p11b

ggsave("../figs/hlm_movement_variability.png",width = 7 ,height=7,dpi = 300)

p11b_M10 <- ggplot(ee2%>%filter(var_type == "volatility", movement_type == "M10"), aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank(),
        axis.text.x = element_text(angle=45))+
  labs(title = "",y="daytime movement volatility",x="")+
  scale_y_continuous(limits = c(-5.8,-3.8),breaks = seq(-5.8,-3.8, by = 0.4))+
  stat_pvalue_manual(stat.test%>%filter(var_type == "volatility", movement_type == "M10"),label = "p_sig",vjust = 1)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))

time

raw value

model_formula_mean <- val ~ Allocation * time_type * phase + posLog + pos_vmu + Age + Gender + (time_type|ID)

df_hlm_lmer_time_mean <- df_hlm%>%
  select(ID,Allocation,phase,day,time_L5_mean,time_M10_mean,Age,Gender,posLog,negLog, pos_vmu)%>%
  pivot_longer(c(time_L5_mean,time_M10_mean),names_pattern =  "time_(.*)_mean",names_to = c("time_type"),values_to = "val")

df_hlm_lmer_time_mean$time_type <- factor(df_hlm_lmer_time_mean$time_type)
df_hlm_lmer_time_mean$time_type <- relevel(df_hlm_lmer_time_mean$time_type, ref = "M10")

mod_time_mean_M10 <- lmer(model_formula_mean,data = df_hlm_lmer_time_mean)
kable(round(summary(mod_time_mean_M10)$coefficients,3))

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	14.566	1.340	100.110	10.869	0.000
Allocationlithium	0.613	0.746	52.525	0.822	0.415
time_typeL5	13.661	0.451	192.100	30.298	0.000
phasepost1	-0.756	0.460	1086.546	-1.645	0.100
phasepost2	-0.468	0.494	1099.553	-0.946	0.344
phasepost3	-0.018	0.455	1098.524	-0.039	0.969
phasepost4	1.022	0.505	1095.268	2.022	0.043
posLog	-0.477	0.307	989.204	-1.557	0.120
pos_vmu	0.036	0.093	996.881	0.384	0.701
Age	-0.038	0.024	29.776	-1.599	0.120
GenderM	0.058	0.567	31.615	0.103	0.919
Allocationlithium:time_typeL5	-1.164	0.637	218.428	-1.829	0.069
Allocationlithium:phasepost1	-0.504	0.675	1095.374	-0.747	0.455
Allocationlithium:phasepost2	-1.230	0.682	1097.438	-1.803	0.072
Allocationlithium:phasepost3	-1.547	0.658	1098.210	-2.350	0.019
Allocationlithium:phasepost4	-1.691	0.713	1097.265	-2.370	0.018
time_typeL5:phasepost1	0.341	0.642	1091.042	0.532	0.595
time_typeL5:phasepost2	-0.244	0.683	1084.170	-0.358	0.721
time_typeL5:phasepost3	-0.612	0.629	1076.765	-0.974	0.330
time_typeL5:phasepost4	-1.503	0.690	1021.463	-2.177	0.030
Allocationlithium:time_typeL5:phasepost1	0.949	0.943	1093.503	1.007	0.314
Allocationlithium:time_typeL5:phasepost2	1.569	0.948	1083.099	1.655	0.098
Allocationlithium:time_typeL5:phasepost3	1.918	0.907	1059.946	2.116	0.035
Allocationlithium:time_typeL5:phasepost4	1.520	0.973	1041.081	1.562	0.119

df_hlm_lmer_time_mean$time_type <- factor(df_hlm_lmer_time_mean$time_type)
df_hlm_lmer_time_mean$time_type <- relevel(df_hlm_lmer_time_mean$time_type, ref = "L5")

mod_time_mean_L5 <- lmer(model_formula_mean,data = df_hlm_lmer_time_mean)
kable(round(summary(mod_time_mean_L5)$coefficients,3))

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	28.228	1.321	96.001	21.366	0.000
Allocationlithium	-0.551	0.678	58.541	-0.812	0.420
time_typeM10	-13.661	0.451	192.105	-30.298	0.000
phasepost1	-0.414	0.459	1088.401	-0.903	0.367
phasepost2	-0.712	0.492	1099.492	-1.446	0.148
phasepost3	-0.630	0.453	1098.461	-1.390	0.165
phasepost4	-0.482	0.501	1079.869	-0.961	0.337
posLog	-0.477	0.307	989.206	-1.557	0.120
pos_vmu	0.036	0.093	996.883	0.384	0.701
Age	-0.038	0.024	29.777	-1.599	0.120
GenderM	0.058	0.567	31.616	0.103	0.919
Allocationlithium:time_typeM10	1.164	0.637	218.433	1.829	0.069
Allocationlithium:phasepost1	0.445	0.673	1097.623	0.662	0.508
Allocationlithium:phasepost2	0.339	0.680	1098.723	0.499	0.618
Allocationlithium:phasepost3	0.371	0.655	1096.811	0.567	0.571
Allocationlithium:phasepost4	-0.171	0.708	1084.489	-0.241	0.809
time_typeM10:phasepost1	-0.341	0.642	1091.042	-0.532	0.595
time_typeM10:phasepost2	0.244	0.683	1084.170	0.358	0.721
time_typeM10:phasepost3	0.612	0.629	1076.764	0.974	0.330
time_typeM10:phasepost4	1.503	0.690	1021.464	2.177	0.030
Allocationlithium:time_typeM10:phasepost1	-0.949	0.943	1093.503	-1.007	0.314
Allocationlithium:time_typeM10:phasepost2	-1.569	0.948	1083.098	-1.655	0.098
Allocationlithium:time_typeM10:phasepost3	-1.918	0.907	1059.946	-2.116	0.035
Allocationlithium:time_typeM10:phasepost4	-1.520	0.973	1041.081	-1.562	0.119

plot_model(mod_time_mean_M10, type = "pred", terms = c("phase", "Allocation","time_type"))+
  scale_color_manual(values=orangeNavyPalette)+
  theme_bw()+
  labs(title = "predicted value of onset time",y="onset time")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

df_hlm_lmer_time_mean$time_type <- relevel(df_hlm_lmer_time_mean$time_type, ref = "M10")
ee3 <- as.data.frame(Effect(c("phase", "Allocation","time_type"),mod_time_mean_M10))

stat.test_M10 <- as.data.frame(round(summary(mod_time_mean_M10)$coefficients,3))[13:16,]
stat.test_L5 <- as.data.frame(round(summary(mod_time_mean_L5)$coefficients,3))[13:16,]
stat.test <- rbind(stat.test_L5,stat.test_M10)%>%rename(p_val = 'Pr(>|t|)')
stat.test[1:4,"time_type"] <- "L5" 
stat.test[5:8,"time_type"] <- "M10" 

stat.test$group1 <- c("post1","post2","post3","post4")
stat.test$group2 <- c("post1","post2","post3","post4")

stat.test <- left_join(stat.test,ee3[ee3$phase != "pre" & ee3$Allocation == "control",]%>%
  rename(group1 = phase))%>%
  mutate(y.position = if_else(time_type == "L5",ee3$upper[11],if_else(time_type == "M10", ee3$upper[1],0)))%>%
  mutate(p_sig = if_else(p_val < 0.001,"***", if_else(p_val < 0.01, "**", if_else(p_val < 0.05, "*", "ns")) ))

## Joining, by = c("time_type", "group1")

p4b <- ggplot(ee3, aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  facet_grid(~time_type)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank(),
        axis.text.x = element_text(angle=45))+
  labs(title = "",y="movement onset time (hrs)")+ 
  scale_y_continuous(limits = c(10,30),breaks = seq(10,30, by = 2.5))+
  stat_pvalue_manual(stat.test, label = "p_sig",vjust = -2)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))
p4b

ggsave("../figs/hlm_time_val.png",width = 6 ,height=4,dpi = 300)

# p4a/p4b
# ggsave("../figs/fig2_vals.png",width = 7 ,height=7,dpi = 300)

p4b_M10 <- ggplot(ee3%>%filter(time_type == "M10"), aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank(),
        axis.text.x = element_text(angle=45))+
  labs(title = "",y="daytime movement onset time",x="")+ 
  scale_y_continuous(limits = c(8,16),breaks = seq(8, 16, by = 2))+
  stat_pvalue_manual(stat.test%>%filter(time_type == "M10"), label = "p_sig",vjust = -4)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))
p4b_M10

variability

model_formula <- val ~ Allocation * var_type * time_type * phase + posLog + pos_vmu + Age + Gender + (var_type * time_type|ID)

df_hlm_lmer_time_variability <- df_hlm%>%
  select(ID,Allocation,day,phase,time_L5_s,time_L5_vmu, time_M10_s,time_M10_vmu,Age,Gender,posLog, negLog, pos_vmu)%>%
  pivot_longer(c(time_L5_s,time_L5_vmu, time_M10_s,time_M10_vmu),names_pattern =  "time_(.*)_(.*)",names_to = c("time_type","var_type"),values_to = "val")%>%
  mutate(var_type = if_else(var_type == "vmu","volatility","noise"))

df_hlm_lmer_time_variability$var_type <- factor(df_hlm_lmer_time_variability$var_type)
df_hlm_lmer_time_variability$var_type <- relevel(df_hlm_lmer_time_variability$var_type, ref = "noise")
df_hlm_lmer_time_variability$time_type <- factor(df_hlm_lmer_time_variability$time_type)
df_hlm_lmer_time_variability$time_type <- relevel(df_hlm_lmer_time_variability$time_type, ref = "M10")

mod_time_var_s_M10 <- lmer(model_formula,data = df_hlm_lmer_time_variability)
kable(round(summary(mod_time_var_s_M10)$coefficients,3))

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	-1.274	0.128	95.836	-9.954	0.000
Allocationlithium	0.128	0.093	58.127	1.377	0.174
var_typevolatility	-2.649	0.181	37.026	-14.605	0.000
time_typeL5	-0.155	0.103	52.298	-1.497	0.140
phasepost1	0.044	0.058	3382.034	0.747	0.455
phasepost2	0.052	0.060	3397.474	0.864	0.387
phasepost3	0.122	0.060	3393.539	2.047	0.041
phasepost4	0.138	0.061	3392.187	2.268	0.023
posLog	0.068	0.025	2699.760	2.704	0.007
pos_vmu	0.002	0.008	2335.667	0.198	0.843
Age	-0.005	0.002	28.986	-2.058	0.049
GenderM	0.019	0.054	31.386	0.358	0.723
Allocationlithium:var_typevolatility	-0.203	0.244	37.643	-0.834	0.409
Allocationlithium:time_typeL5	-0.110	0.140	54.807	-0.787	0.435
var_typevolatility:time_typeL5	-0.115	0.174	44.472	-0.660	0.513
Allocationlithium:phasepost1	-0.035	0.080	3384.337	-0.442	0.659
Allocationlithium:phasepost2	-0.092	0.082	3388.245	-1.118	0.264
Allocationlithium:phasepost3	-0.113	0.084	3392.596	-1.343	0.179
Allocationlithium:phasepost4	-0.048	0.085	3396.764	-0.565	0.572
var_typevolatility:phasepost1	-0.661	0.082	3373.348	-8.024	0.000
var_typevolatility:phasepost2	-0.962	0.084	3383.992	-11.407	0.000
var_typevolatility:phasepost3	-1.398	0.084	3381.797	-16.621	0.000
var_typevolatility:phasepost4	-1.557	0.086	3381.627	-18.142	0.000
time_typeL5:phasepost1	-0.251	0.082	3376.343	-3.045	0.002
time_typeL5:phasepost2	-0.237	0.084	3392.179	-2.811	0.005
time_typeL5:phasepost3	-0.367	0.084	3389.237	-4.369	0.000
time_typeL5:phasepost4	-0.408	0.086	3389.084	-4.754	0.000
Allocationlithium:var_typevolatility:time_typeL5	0.186	0.235	46.035	0.792	0.432
Allocationlithium:var_typevolatility:phasepost1	0.171	0.113	3377.229	1.514	0.130
Allocationlithium:var_typevolatility:phasepost2	0.237	0.117	3379.539	2.030	0.042
Allocationlithium:var_typevolatility:phasepost3	0.458	0.119	3381.999	3.858	0.000
Allocationlithium:var_typevolatility:phasepost4	0.381	0.120	3385.310	3.164	0.002
Allocationlithium:time_typeL5:phasepost1	0.036	0.113	3382.715	0.320	0.749
Allocationlithium:time_typeL5:phasepost2	-0.028	0.117	3386.073	-0.244	0.807
Allocationlithium:time_typeL5:phasepost3	0.061	0.118	3389.866	0.512	0.608
Allocationlithium:time_typeL5:phasepost4	0.009	0.120	3394.365	0.073	0.942
var_typevolatility:time_typeL5:phasepost1	0.032	0.116	3374.340	0.274	0.784
var_typevolatility:time_typeL5:phasepost2	0.110	0.119	3387.403	0.926	0.354
var_typevolatility:time_typeL5:phasepost3	0.696	0.119	3384.923	5.857	0.000
var_typevolatility:time_typeL5:phasepost4	0.914	0.121	3384.781	7.529	0.000
Allocationlithium:var_typevolatility:time_typeL5:phasepost1	0.033	0.160	3379.339	0.204	0.838
Allocationlithium:var_typevolatility:time_typeL5:phasepost2	0.257	0.165	3382.111	1.558	0.119
Allocationlithium:var_typevolatility:time_typeL5:phasepost3	-0.377	0.168	3385.311	-2.247	0.025
Allocationlithium:var_typevolatility:time_typeL5:phasepost4	-0.269	0.170	3389.410	-1.579	0.114

df_hlm_lmer_time_variability$var_type <- relevel(df_hlm_lmer_time_variability$var_type, ref = "volatility")
df_hlm_lmer_time_variability$time_type <- relevel(df_hlm_lmer_time_variability$time_type, ref = "M10")

#same model, different reference level
mod_time_var_vmu_M10 <- lmer(model_formula,data = df_hlm_lmer_time_variability)
kable(round(summary(mod_time_var_vmu_M10)$coefficients,3))

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	-3.923	0.187	69.874	-21.007	0.000
Allocationlithium	-0.076	0.204	36.120	-0.371	0.713
var_typenoise	2.649	0.181	37.057	14.609	0.000
time_typeL5	-0.270	0.138	41.058	-1.954	0.058
phasepost1	-0.617	0.058	3376.304	-10.553	0.000
phasepost2	-0.910	0.060	3383.686	-15.168	0.000
phasepost3	-1.276	0.060	3380.657	-21.348	0.000
phasepost4	-1.420	0.061	3379.188	-23.302	0.000
posLog	0.068	0.025	2699.634	2.704	0.007
pos_vmu	0.002	0.008	2335.506	0.198	0.843
Age	-0.005	0.002	28.990	-2.058	0.049
GenderM	0.019	0.054	31.390	0.357	0.723
Allocationlithium:var_typenoise	0.203	0.244	37.674	0.834	0.409
Allocationlithium:time_typeL5	0.076	0.186	42.229	0.408	0.685
var_typenoise:time_typeL5	0.115	0.174	44.494	0.660	0.513
Allocationlithium:phasepost1	0.136	0.080	3373.991	1.699	0.089
Allocationlithium:phasepost2	0.145	0.083	3375.669	1.751	0.080
Allocationlithium:phasepost3	0.345	0.084	3378.680	4.099	0.000
Allocationlithium:phasepost4	0.333	0.086	3381.195	3.893	0.000
var_typenoise:phasepost1	0.661	0.082	3373.330	8.024	0.000
var_typenoise:phasepost2	0.962	0.084	3383.974	11.407	0.000
var_typenoise:phasepost3	1.398	0.084	3381.779	16.621	0.000
var_typenoise:phasepost4	1.557	0.086	3381.610	18.141	0.000
time_typeL5:phasepost1	-0.219	0.082	3372.806	-2.657	0.008
time_typeL5:phasepost2	-0.126	0.084	3383.055	-1.498	0.134
time_typeL5:phasepost3	0.329	0.084	3381.127	3.914	0.000
time_typeL5:phasepost4	0.506	0.086	3380.988	5.893	0.000
Allocationlithium:var_typenoise:time_typeL5	-0.186	0.235	46.058	-0.792	0.432
Allocationlithium:var_typenoise:phasepost1	-0.171	0.113	3377.211	-1.514	0.130
Allocationlithium:var_typenoise:phasepost2	-0.237	0.117	3379.521	-2.030	0.042
Allocationlithium:var_typenoise:phasepost3	-0.458	0.119	3381.982	-3.858	0.000
Allocationlithium:var_typenoise:phasepost4	-0.381	0.120	3385.293	-3.164	0.002
Allocationlithium:time_typeL5:phasepost1	0.069	0.113	3376.569	0.608	0.543
Allocationlithium:time_typeL5:phasepost2	0.228	0.117	3378.776	1.959	0.050
Allocationlithium:time_typeL5:phasepost3	-0.316	0.119	3381.270	-2.664	0.008
Allocationlithium:time_typeL5:phasepost4	-0.260	0.120	3384.668	-2.159	0.031
var_typenoise:time_typeL5:phasepost1	-0.032	0.116	3374.322	-0.274	0.784
var_typenoise:time_typeL5:phasepost2	-0.110	0.119	3387.385	-0.926	0.354
var_typenoise:time_typeL5:phasepost3	-0.696	0.119	3384.905	-5.857	0.000
var_typenoise:time_typeL5:phasepost4	-0.914	0.121	3384.763	-7.529	0.000
Allocationlithium:var_typenoise:time_typeL5:phasepost1	-0.033	0.160	3379.321	-0.204	0.838
Allocationlithium:var_typenoise:time_typeL5:phasepost2	-0.257	0.165	3382.093	-1.558	0.119
Allocationlithium:var_typenoise:time_typeL5:phasepost3	0.377	0.168	3385.293	2.247	0.025
Allocationlithium:var_typenoise:time_typeL5:phasepost4	0.269	0.170	3389.392	1.579	0.114

df_hlm_lmer_time_variability$Allocation <- factor(df_hlm_lmer_time_variability$Allocation, levels = c("lithium","control"))
kable(round(summary(lmer(model_formula,data = df_hlm_lmer_time_variability))$coefficients,3))

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	-3.999	0.169	76.963	-23.623	0.000
Allocationcontrol	0.076	0.204	36.098	0.370	0.713
var_typenoise	2.852	0.163	38.441	17.534	0.000
time_typeL5	-0.194	0.125	43.712	-1.554	0.127
phasepost1	-0.481	0.055	3379.372	-8.729	0.000
phasepost2	-0.766	0.057	3376.181	-13.391	0.000
phasepost3	-0.931	0.059	3377.152	-15.724	0.000
phasepost4	-1.087	0.060	3382.310	-18.151	0.000
posLog	0.068	0.025	2699.721	2.704	0.007
pos_vmu	0.002	0.008	2335.612	0.198	0.843
Age	-0.005	0.002	28.989	-2.058	0.049
GenderM	0.019	0.054	31.389	0.358	0.723
Allocationcontrol:var_typenoise	-0.203	0.244	37.652	-0.834	0.409
Allocationcontrol:time_typeL5	-0.076	0.186	42.207	-0.408	0.685
var_typenoise:time_typeL5	-0.071	0.158	48.058	-0.452	0.653
Allocationcontrol:phasepost1	-0.136	0.080	3374.004	-1.699	0.089
Allocationcontrol:phasepost2	-0.145	0.083	3375.682	-1.751	0.080
Allocationcontrol:phasepost3	-0.345	0.084	3378.693	-4.099	0.000
Allocationcontrol:phasepost4	-0.333	0.086	3381.208	-3.893	0.000
var_typenoise:phasepost1	0.490	0.078	3381.263	6.303	0.000
var_typenoise:phasepost2	0.725	0.081	3374.115	9.010	0.000
var_typenoise:phasepost3	0.940	0.084	3382.479	11.241	0.000
var_typenoise:phasepost4	1.176	0.084	3389.126	13.924	0.000
time_typeL5:phasepost1	-0.150	0.078	3380.515	-1.932	0.053
time_typeL5:phasepost2	0.102	0.081	3373.612	1.267	0.205
time_typeL5:phasepost3	0.013	0.084	3381.464	0.158	0.874
time_typeL5:phasepost4	0.246	0.085	3388.174	2.910	0.004
Allocationcontrol:var_typenoise:time_typeL5	0.186	0.235	46.042	0.792	0.432
Allocationcontrol:var_typenoise:phasepost1	0.171	0.113	3377.226	1.514	0.130
Allocationcontrol:var_typenoise:phasepost2	0.237	0.117	3379.535	2.030	0.042
Allocationcontrol:var_typenoise:phasepost3	0.458	0.119	3381.996	3.858	0.000
Allocationcontrol:var_typenoise:phasepost4	0.381	0.120	3385.307	3.164	0.002
Allocationcontrol:time_typeL5:phasepost1	-0.069	0.113	3376.582	-0.608	0.543
Allocationcontrol:time_typeL5:phasepost2	-0.228	0.117	3378.789	-1.959	0.050
Allocationcontrol:time_typeL5:phasepost3	0.316	0.119	3381.282	2.664	0.008
Allocationcontrol:time_typeL5:phasepost4	0.260	0.120	3384.680	2.159	0.031
var_typenoise:time_typeL5:phasepost1	-0.065	0.110	3384.360	-0.588	0.557
var_typenoise:time_typeL5:phasepost2	-0.367	0.114	3375.348	-3.227	0.001
var_typenoise:time_typeL5:phasepost3	-0.320	0.118	3385.693	-2.703	0.007
var_typenoise:time_typeL5:phasepost4	-0.645	0.119	3393.402	-5.401	0.000
Allocationcontrol:var_typenoise:time_typeL5:phasepost1	0.033	0.160	3379.334	0.204	0.838
Allocationcontrol:var_typenoise:time_typeL5:phasepost2	0.257	0.165	3382.107	1.558	0.119
Allocationcontrol:var_typenoise:time_typeL5:phasepost3	-0.377	0.168	3385.306	-2.247	0.025
Allocationcontrol:var_typenoise:time_typeL5:phasepost4	-0.269	0.170	3389.406	-1.579	0.114

df_hlm_lmer_time_variability$Allocation <- relevel(df_hlm_lmer_time_variability$Allocation, ref = "control")

# #same model, different reference level
df_hlm_lmer_time_variability$var_type <- relevel(df_hlm_lmer_time_variability$var_type, ref = "noise")
df_hlm_lmer_time_variability$time_type <- relevel(df_hlm_lmer_time_variability$time_type, ref = "L5")

mod_time_var_s_L5 <- lmer(model_formula,data = df_hlm_lmer_time_variability)
kable(round(summary(mod_time_var_s_L5)$coefficients,3))

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	-1.429	0.133	99.677	-10.773	0.000
Allocationlithium	0.017	0.104	50.012	0.167	0.868
var_typevolatility	-2.764	0.146	40.036	-18.941	0.000
time_typeM10	0.155	0.103	52.278	1.497	0.140
phasepost1	-0.207	0.058	3380.288	-3.542	0.000
phasepost2	-0.185	0.060	3395.223	-3.089	0.002
phasepost3	-0.245	0.060	3390.979	-4.104	0.000
phasepost4	-0.270	0.061	3389.483	-4.435	0.000
posLog	0.068	0.025	2699.736	2.703	0.007
pos_vmu	0.002	0.008	2335.626	0.198	0.843
Age	-0.005	0.002	28.991	-2.058	0.049
GenderM	0.019	0.054	31.391	0.358	0.723
Allocationlithium:var_typevolatility	-0.017	0.196	41.048	-0.087	0.931
Allocationlithium:time_typeM10	0.110	0.140	54.786	0.787	0.435
var_typevolatility:time_typeM10	0.115	0.174	44.465	0.659	0.513
Allocationlithium:phasepost1	0.001	0.080	3381.506	0.011	0.992
Allocationlithium:phasepost2	-0.121	0.082	3384.946	-1.463	0.144
Allocationlithium:phasepost3	-0.052	0.084	3389.551	-0.619	0.536
Allocationlithium:phasepost4	-0.039	0.085	3393.976	-0.461	0.645
var_typevolatility:phasepost1	-0.629	0.082	3374.633	-7.637	0.000
var_typevolatility:phasepost2	-0.852	0.084	3388.044	-10.101	0.000
var_typevolatility:phasepost3	-0.701	0.084	3385.452	-8.344	0.000
var_typevolatility:phasepost4	-0.644	0.086	3385.276	-7.500	0.000
time_typeM10:phasepost1	0.251	0.082	3376.346	3.045	0.002
time_typeM10:phasepost2	0.237	0.084	3392.181	2.811	0.005
time_typeM10:phasepost3	0.367	0.084	3389.240	4.369	0.000
time_typeM10:phasepost4	0.408	0.086	3389.087	4.755	0.000
Allocationlithium:var_typevolatility:time_typeM10	-0.186	0.235	46.028	-0.792	0.432
Allocationlithium:var_typevolatility:phasepost1	0.204	0.113	3379.748	1.803	0.071
Allocationlithium:var_typevolatility:phasepost2	0.494	0.117	3382.631	4.234	0.000
Allocationlithium:var_typevolatility:phasepost3	0.081	0.119	3385.807	0.682	0.495
Allocationlithium:var_typevolatility:phasepost4	0.112	0.120	3389.898	0.931	0.352
Allocationlithium:time_typeM10:phasepost1	-0.036	0.113	3382.718	-0.320	0.749
Allocationlithium:time_typeM10:phasepost2	0.028	0.117	3386.075	0.244	0.807
Allocationlithium:time_typeM10:phasepost3	-0.061	0.118	3389.869	-0.512	0.608
Allocationlithium:time_typeM10:phasepost4	-0.009	0.120	3394.367	-0.073	0.942
var_typevolatility:time_typeM10:phasepost1	-0.032	0.116	3374.343	-0.274	0.784
var_typevolatility:time_typeM10:phasepost2	-0.110	0.119	3387.406	-0.926	0.354
var_typevolatility:time_typeM10:phasepost3	-0.696	0.119	3384.926	-5.857	0.000
var_typevolatility:time_typeM10:phasepost4	-0.914	0.121	3384.784	-7.529	0.000
Allocationlithium:var_typevolatility:time_typeM10:phasepost1	-0.033	0.160	3379.342	-0.204	0.838
Allocationlithium:var_typevolatility:time_typeM10:phasepost2	-0.257	0.165	3382.114	-1.558	0.119
Allocationlithium:var_typevolatility:time_typeM10:phasepost3	0.377	0.168	3385.314	2.247	0.025
Allocationlithium:var_typevolatility:time_typeM10:phasepost4	0.269	0.170	3389.413	1.579	0.114

df_hlm_lmer_time_variability$var_type <- relevel(df_hlm_lmer_time_variability$var_type, ref = "volatility")

df_hlm_lmer_time_variability$time_type <- relevel(df_hlm_lmer_time_variability$time_type, ref = "L5")
mod_time_var_vmu_L5 <- lmer(model_formula,data = df_hlm_lmer_time_variability)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00215358 (tol = 0.002, component 1)

kable(round(summary(mod_time_var_vmu_L5)$coefficients,3))

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	-4.193	0.183	71.618	-22.874	0.000
Allocationlithium	0.000	0.199	36.374	0.001	0.999
var_typenoise	2.764	0.146	40.064	18.945	0.000
time_typeM10	0.270	0.138	41.047	1.953	0.058
phasepost1	-0.836	0.058	3376.356	-14.295	0.000
phasepost2	-1.037	0.060	3384.210	-17.275	0.000
phasepost3	-0.946	0.060	3381.137	-15.839	0.000
phasepost4	-0.914	0.061	3379.763	-14.997	0.000
posLog	0.068	0.025	2699.711	2.703	0.007
pos_vmu	0.002	0.008	2335.602	0.198	0.843
Age	-0.005	0.002	28.991	-2.058	0.049
GenderM	0.019	0.054	31.391	0.358	0.723
Allocationlithium:var_typenoise	0.017	0.196	41.077	0.087	0.931
Allocationlithium:time_typeM10	-0.076	0.186	42.218	-0.408	0.685
var_typenoise:time_typeM10	-0.115	0.174	44.490	-0.660	0.513
Allocationlithium:phasepost1	0.205	0.080	3374.255	2.559	0.011
Allocationlithium:phasepost2	0.373	0.083	3375.930	4.519	0.000
Allocationlithium:phasepost3	0.029	0.084	3379.156	0.343	0.732
Allocationlithium:phasepost4	0.073	0.085	3381.853	0.851	0.395
var_typenoise:phasepost1	0.629	0.082	3374.619	7.637	0.000
var_typenoise:phasepost2	0.852	0.084	3388.031	10.101	0.000
var_typenoise:phasepost3	0.701	0.084	3385.439	8.344	0.000
var_typenoise:phasepost4	0.644	0.086	3385.263	7.500	0.000
time_typeM10:phasepost1	0.219	0.082	3372.812	2.657	0.008
time_typeM10:phasepost2	0.126	0.084	3383.060	1.498	0.134
time_typeM10:phasepost3	-0.329	0.084	3381.133	-3.914	0.000
time_typeM10:phasepost4	-0.506	0.086	3380.994	-5.893	0.000
Allocationlithium:var_typenoise:time_typeM10	0.186	0.235	46.054	0.792	0.432
Allocationlithium:var_typenoise:phasepost1	-0.204	0.113	3379.735	-1.803	0.071
Allocationlithium:var_typenoise:phasepost2	-0.494	0.117	3382.618	-4.234	0.000
Allocationlithium:var_typenoise:phasepost3	-0.081	0.119	3385.794	-0.682	0.495
Allocationlithium:var_typenoise:phasepost4	-0.112	0.120	3389.886	-0.931	0.352
Allocationlithium:time_typeM10:phasepost1	-0.069	0.113	3376.574	-0.608	0.543
Allocationlithium:time_typeM10:phasepost2	-0.228	0.117	3378.782	-1.959	0.050
Allocationlithium:time_typeM10:phasepost3	0.316	0.119	3381.275	2.664	0.008
Allocationlithium:time_typeM10:phasepost4	0.260	0.120	3384.673	2.159	0.031
var_typenoise:time_typeM10:phasepost1	0.032	0.116	3374.329	0.274	0.784
var_typenoise:time_typeM10:phasepost2	0.110	0.119	3387.394	0.926	0.354
var_typenoise:time_typeM10:phasepost3	0.696	0.119	3384.914	5.857	0.000
var_typenoise:time_typeM10:phasepost4	0.914	0.121	3384.771	7.529	0.000
Allocationlithium:var_typenoise:time_typeM10:phasepost1	0.033	0.160	3379.328	0.204	0.838
Allocationlithium:var_typenoise:time_typeM10:phasepost2	0.257	0.165	3382.101	1.558	0.119
Allocationlithium:var_typenoise:time_typeM10:phasepost3	-0.377	0.168	3385.301	-2.247	0.025
Allocationlithium:var_typenoise:time_typeM10:phasepost4	-0.269	0.170	3389.402	-1.579	0.114

plot_model(mod_time_var_s_M10, type = "pred", terms = c("phase", "Allocation","time_type","var_type"))+
  labs(title="predicted value of onset time variability",y="noise")

ggsave("../figs/hlm_time_variability_comparison.png",width = 7 ,height=7,dpi = 300)

df_hlm_lmer_time_variability$var_type <- relevel(df_hlm_lmer_time_variability$var_type, ref = "volatility")
df_hlm_lmer_time_variability$time_type <- relevel(df_hlm_lmer_time_variability$time_type, ref = "M10")
ee4 <- as.data.frame(Effect(c("phase", "Allocation","time_type","var_type"),mod_time_var_vmu_M10))

stat.test_s_M10 <- as.data.frame(round(summary(mod_time_var_s_M10)$coefficients,3))[16:19,]
stat.test_s_L5 <- as.data.frame(round(summary(mod_time_var_s_L5)$coefficients,3))[16:19,]
stat.test_vmu_M10 <- as.data.frame(round(summary(mod_time_var_vmu_M10)$coefficients,3))[16:19,]
stat.test_vmu_L5 <- as.data.frame(round(summary(mod_time_var_vmu_L5)$coefficients,3))[16:19,]

stat.test <- Reduce(rbind,list(stat.test_s_M10,stat.test_s_L5,stat.test_vmu_M10,stat.test_vmu_L5))%>%rename(p_val = 'Pr(>|t|)')
stat.test[c(1:4,9:12),"time_type"] <- "M10" 
stat.test[c(5:8,13:16),"time_type"] <- "L5" 

stat.test[1:8,"var_type"] <- "noise" 
stat.test[9:16,"var_type"] <- "volatility" 

stat.test$group1 <- c("post1","post2","post3","post4")
stat.test$group2 <- c("post1","post2","post3","post4")

stat.test <- left_join(stat.test,ee4[ee4$phase != "pre" & ee4$Allocation == "lithium",]%>%
  rename(group1 = phase), by = c("time_type","var_type","group1"))%>%
  mutate(y.position = if_else(var_type == "volatility" & time_type == "L5",ee4$upper[11],if_else(var_type == "volatility" & time_type == "M10", ee4$upper[1],if_else(var_type == "noise" & time_type == "L5",ee4$upper[36],if_else(var_type == "noise" & time_type == "M10",ee4$upper[21],0)))))%>%
  mutate(p_sig = if_else(p_val < 0.001,"***", if_else(p_val < 0.01, "**", if_else(p_val < 0.05, "*", "ns")) ))

#plot noise first
p12a <- ggplot(ee4%>%filter(var_type == "noise"), aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  facet_wrap(~time_type)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank())+
  labs(title = "",y="movement noise")+
  stat_pvalue_manual(stat.test%>%filter(var_type == "noise"),label = "p_sig",vjust = -1)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))

#plot volatility first
p12b <- ggplot(ee4%>%filter(var_type == "volatility"), aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  facet_wrap(~time_type)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank())+
  labs(title = "",y="movement volatility")+
  stat_pvalue_manual(stat.test%>%filter(var_type == "volatility"),label = "p_sig",vjust = 1)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))

p12a/p12b

ggsave("../figs/hlm_time_variability.png",width = 7 ,height=7,dpi = 300)

p12b_M10 <- ggplot(ee4%>%filter(var_type == "volatility", time_type == "M10"), aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank(),
        axis.text.x = element_text(angle=45))+
  labs(title = "",y="daytime onset time volatility",x="")+
  stat_pvalue_manual(stat.test%>%filter(var_type == "volatility", time_type == "M10"),label = "p_sig",vjust = 1)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))

#manuscript graphs

p4a_M10 + p11b_M10+
  plot_annotation(tag_levels = c("A"))

ggsave("../figs/movement.png",width = 7 ,height=3,dpi = 300)

p4b_M10 + p12b_M10+
  plot_annotation(tag_levels = c("A"))

ggsave("../figs/time.png",width = 7 ,height=3,dpi = 300)

(p4a_M10 + p11b_M10) / (p4b_M10 + p12b_M10) +
  plot_annotation(tag_levels = c("A"))

ggsave("../figs/FIG4_movement_n_time.png",width = 7 ,height=6,dpi = 300)
  
(p8 + p9) / (p4a_M10 + p11b_M10) / (p4b_M10 + p12b_M10) +
  plot_annotation(tag_levels = c("A"))

ggsave("../figs/movement_n_time_n_pa.png",width = 7 ,height=10,dpi = 300)

correlating change in vmu and mean

df_hlm$phase <- factor(df_hlm$phase, labels = c("run-in","week 1","week 2","week 3","week 4"))

cA1 <- lmer(time_M10_mean ~ Allocation * movement_M10_mean * phase + Age + Gender + (1|ID),data = df_hlm)
summary(cA1)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: time_M10_mean ~ Allocation * movement_M10_mean * phase + Age +  
##     Gender + (1 | ID)
##    Data: df_hlm
## 
## REML criterion at convergence: 4133.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3191 -0.7139 -0.0724  0.7163  2.8656 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept)  2.349   1.533   
##  Residual             11.257   3.355   
## Number of obs: 767, groups:  ID, 34
## 
## Fixed effects:
##                                                   Estimate Std. Error
## (Intercept)                                      14.083237   1.491655
## Allocationlithium                                -0.306840   1.432668
## movement_M10_mean                                -0.014464   0.017957
## phaseweek 1                                      -1.711169   1.547788
## phaseweek 2                                      -2.010896   1.589343
## phaseweek 3                                      -2.235538   1.399565
## phaseweek 4                                      -1.483085   1.616221
## Age                                              -0.040895   0.025603
## GenderM                                           0.803187   0.603491
## Allocationlithium:movement_M10_mean               0.006565   0.021387
## Allocationlithium:phaseweek 1                     1.473474   2.040221
## Allocationlithium:phaseweek 2                     1.561965   1.925902
## Allocationlithium:phaseweek 3                     1.316191   1.811787
## Allocationlithium:phaseweek 4                     0.924409   2.045464
## movement_M10_mean:phaseweek 1                     0.012929   0.025778
## movement_M10_mean:phaseweek 2                     0.026941   0.024114
## movement_M10_mean:phaseweek 3                     0.031050   0.021608
## movement_M10_mean:phaseweek 4                     0.038700   0.025195
## Allocationlithium:movement_M10_mean:phaseweek 1  -0.028480   0.033281
## Allocationlithium:movement_M10_mean:phaseweek 2  -0.040369   0.030056
## Allocationlithium:movement_M10_mean:phaseweek 3  -0.035999   0.029345
## Allocationlithium:movement_M10_mean:phaseweek 4  -0.035553   0.033116
##                                                         df t value Pr(>|t|)    
## (Intercept)                                     110.389365   9.441 7.25e-16 ***
## Allocationlithium                               387.029151  -0.214    0.831    
## movement_M10_mean                               723.112305  -0.805    0.421    
## phaseweek 1                                     724.774360  -1.106    0.269    
## phaseweek 2                                     721.807503  -1.265    0.206    
## phaseweek 3                                     732.081078  -1.597    0.111    
## phaseweek 4                                     734.518937  -0.918    0.359    
## Age                                              30.765009  -1.597    0.120    
## GenderM                                          29.895186   1.331    0.193    
## Allocationlithium:movement_M10_mean             738.006859   0.307    0.759    
## Allocationlithium:phaseweek 1                   729.009206   0.722    0.470    
## Allocationlithium:phaseweek 2                   722.846139   0.811    0.418    
## Allocationlithium:phaseweek 3                   730.096092   0.726    0.468    
## Allocationlithium:phaseweek 4                   736.339936   0.452    0.651    
## movement_M10_mean:phaseweek 1                   728.339714   0.502    0.616    
## movement_M10_mean:phaseweek 2                   720.840497   1.117    0.264    
## movement_M10_mean:phaseweek 3                   732.617335   1.437    0.151    
## movement_M10_mean:phaseweek 4                   740.436632   1.536    0.125    
## Allocationlithium:movement_M10_mean:phaseweek 1 729.139129  -0.856    0.392    
## Allocationlithium:movement_M10_mean:phaseweek 2 721.536396  -1.343    0.180    
## Allocationlithium:movement_M10_mean:phaseweek 3 729.844951  -1.227    0.220    
## Allocationlithium:movement_M10_mean:phaseweek 4 737.476413  -1.074    0.283    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 22 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

cA2 <- lmer(pos_mean ~ Allocation * movement_M10_mean * phase + Age + Gender + (1|ID),data = df_hlm)
summary(cA2)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pos_mean ~ Allocation * movement_M10_mean * phase + Age + Gender +  
##     (1 | ID)
##    Data: df_hlm
## 
## REML criterion at convergence: 2997
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8470 -0.5705 -0.0628  0.5891  3.2847 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 8.772    2.962   
##  Residual             8.530    2.921   
## Number of obs: 575, groups:  ID, 34
## 
## Fixed effects:
##                                                   Estimate Std. Error
## (Intercept)                                       9.266784   2.098319
## Allocationlithium                                 3.476211   1.653969
## movement_M10_mean                                 0.038176   0.016858
## phaseweek 1                                      -3.554955   1.499697
## phaseweek 2                                      -0.985395   1.493013
## phaseweek 3                                       1.667830   1.305060
## phaseweek 4                                       2.500531   1.477875
## Age                                               0.024768   0.045479
## GenderM                                          -3.611984   1.078351
## Allocationlithium:movement_M10_mean              -0.034361   0.020190
## Allocationlithium:phaseweek 1                    -0.454923   2.004505
## Allocationlithium:phaseweek 2                    -1.433871   1.882547
## Allocationlithium:phaseweek 3                    -3.343205   1.796356
## Allocationlithium:phaseweek 4                    -3.787107   1.995389
## movement_M10_mean:phaseweek 1                     0.034753   0.025166
## movement_M10_mean:phaseweek 2                     0.005019   0.022474
## movement_M10_mean:phaseweek 3                    -0.023569   0.020040
## movement_M10_mean:phaseweek 4                    -0.043633   0.023100
## Allocationlithium:movement_M10_mean:phaseweek 1   0.002307   0.032607
## Allocationlithium:movement_M10_mean:phaseweek 2   0.005677   0.028800
## Allocationlithium:movement_M10_mean:phaseweek 3   0.031761   0.029160
## Allocationlithium:movement_M10_mean:phaseweek 4   0.039144   0.032336
##                                                         df t value Pr(>|t|)    
## (Intercept)                                      50.282666   4.416 5.34e-05 ***
## Allocationlithium                               139.115978   2.102  0.03738 *  
## movement_M10_mean                               544.946598   2.265  0.02393 *  
## phaseweek 1                                     527.556294  -2.370  0.01813 *  
## phaseweek 2                                     526.220935  -0.660  0.50954    
## phaseweek 3                                     529.156148   1.278  0.20182    
## phaseweek 4                                     530.576186   1.692  0.09124 .  
## Age                                              30.619426   0.545  0.58998    
## GenderM                                          30.873126  -3.350  0.00215 ** 
## Allocationlithium:movement_M10_mean             541.477980  -1.702  0.08935 .  
## Allocationlithium:phaseweek 1                   528.689951  -0.227  0.82055    
## Allocationlithium:phaseweek 2                   526.941055  -0.762  0.44660    
## Allocationlithium:phaseweek 3                   528.071771  -1.861  0.06329 .  
## Allocationlithium:phaseweek 4                   534.255909  -1.898  0.05824 .  
## movement_M10_mean:phaseweek 1                   529.065792   1.381  0.16787    
## movement_M10_mean:phaseweek 2                   525.556403   0.223  0.82338    
## movement_M10_mean:phaseweek 3                   528.844678  -1.176  0.24007    
## movement_M10_mean:phaseweek 4                   534.631607  -1.889  0.05945 .  
## Allocationlithium:movement_M10_mean:phaseweek 1 529.057599   0.071  0.94363    
## Allocationlithium:movement_M10_mean:phaseweek 2 525.966282   0.197  0.84382    
## Allocationlithium:movement_M10_mean:phaseweek 3 527.430213   1.089  0.27657    
## Allocationlithium:movement_M10_mean:phaseweek 4 533.493385   1.211  0.22661    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 22 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

cA3 <- lmer(pos_mean ~ Allocation * time_M10_mean * phase + Age + Gender + (1|ID),data = df_hlm)
summary(cA3)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pos_mean ~ Allocation * time_M10_mean * phase + Age + Gender +  
##     (1 | ID)
##    Data: df_hlm
## 
## REML criterion at convergence: 2969.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.63961 -0.57817 -0.07669  0.61344  2.94435 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 8.965    2.994   
##  Residual             8.688    2.948   
## Number of obs: 575, groups:  ID, 34
## 
## Fixed effects:
##                                              Estimate Std. Error        df
## (Intercept)                                  14.34854    2.27998  66.20570
## Allocationlithium                             0.28204    2.07303 263.32115
## time_M10_mean                                -0.22419    0.10113 538.41941
## phaseweek 1                                  -6.73621    1.82388 527.40012
## phaseweek 2                                  -2.87334    1.91234 527.77995
## phaseweek 3                                  -2.31736    1.80492 532.64005
## phaseweek 4                                  -1.98327    1.99983 529.49126
## Age                                           0.01628    0.04601  30.63413
## GenderM                                      -3.54187    1.08844  30.62932
## Allocationlithium:time_M10_mean               0.11098    0.13822 536.02124
## Allocationlithium:phaseweek 1                 2.75382    2.49878 529.07945
## Allocationlithium:phaseweek 2                -0.37921    2.64428 528.12459
## Allocationlithium:phaseweek 3                 2.09386    2.51152 530.81644
## Allocationlithium:phaseweek 4                -2.09369    2.69668 528.82996
## time_M10_mean:phaseweek 1                     0.45386    0.15449 527.13715
## time_M10_mean:phaseweek 2                     0.21242    0.16199 529.56932
## time_M10_mean:phaseweek 3                     0.23268    0.14976 535.76168
## time_M10_mean:phaseweek 4                     0.18096    0.15957 531.87043
## Allocationlithium:time_M10_mean:phaseweek 1  -0.28448    0.20518 528.96994
## Allocationlithium:time_M10_mean:phaseweek 2  -0.09835    0.22198 530.04149
## Allocationlithium:time_M10_mean:phaseweek 3  -0.33724    0.20386 532.29166
## Allocationlithium:time_M10_mean:phaseweek 4   0.02463    0.21174 529.64031
##                                             t value Pr(>|t|)    
## (Intercept)                                   6.293 2.83e-08 ***
## Allocationlithium                             0.136 0.891882    
## time_M10_mean                                -2.217 0.027055 *  
## phaseweek 1                                  -3.693 0.000244 ***
## phaseweek 2                                  -1.503 0.133559    
## phaseweek 3                                  -1.284 0.199730    
## phaseweek 4                                  -0.992 0.321788    
## Age                                           0.354 0.725791    
## GenderM                                      -3.254 0.002773 ** 
## Allocationlithium:time_M10_mean               0.803 0.422361    
## Allocationlithium:phaseweek 1                 1.102 0.270935    
## Allocationlithium:phaseweek 2                -0.143 0.886024    
## Allocationlithium:phaseweek 3                 0.834 0.404825    
## Allocationlithium:phaseweek 4                -0.776 0.437863    
## time_M10_mean:phaseweek 1                     2.938 0.003451 ** 
## time_M10_mean:phaseweek 2                     1.311 0.190330    
## time_M10_mean:phaseweek 3                     1.554 0.120867    
## time_M10_mean:phaseweek 4                     1.134 0.257285    
## Allocationlithium:time_M10_mean:phaseweek 1  -1.387 0.166172    
## Allocationlithium:time_M10_mean:phaseweek 2  -0.443 0.657925    
## Allocationlithium:time_M10_mean:phaseweek 3  -1.654 0.098663 .  
## Allocationlithium:time_M10_mean:phaseweek 4   0.116 0.907456    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 22 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

cB1 <- lmer(time_M10_vmu ~ Allocation * movement_M10_vmu * phase + Age + Gender + (1|ID),data = df_hlm)
summary(cB1)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: time_M10_vmu ~ Allocation * movement_M10_vmu * phase + Age +  
##     Gender + (1 | ID)
##    Data: df_hlm
## 
## REML criterion at convergence: 1700.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1077 -0.6310 -0.0649  0.5670  3.9057 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.3303   0.5747  
##  Residual             0.1951   0.4417  
## Number of obs: 1243, groups:  ID, 34
## 
## Fixed effects:
##                                                  Estimate Std. Error         df
## (Intercept)                                      -2.17590    0.43103   69.79928
## Allocationlithium                                 0.44302    0.37654  286.96038
## movement_M10_vmu                                  0.32468    0.05884 1196.34241
## phaseweek 1                                       0.03898    0.42094 1196.79933
## phaseweek 2                                       0.03499    0.42611 1197.60737
## phaseweek 3                                      -0.13169    0.46659 1198.93349
## phaseweek 4                                      -0.60940    0.48631 1195.77350
## Age                                              -0.01257    0.00858   29.97599
## GenderM                                           0.07032    0.20314   30.10371
## Allocationlithium:movement_M10_vmu                0.14081    0.07449 1199.04329
## Allocationlithium:phaseweek 1                    -1.25961    0.49930 1195.80379
## Allocationlithium:phaseweek 2                    -1.01360    0.51375 1196.67988
## Allocationlithium:phaseweek 3                    -1.76943    0.54808 1198.21475
## Allocationlithium:phaseweek 4                    -1.29420    0.56350 1196.07011
## movement_M10_vmu:phaseweek 1                      0.06931    0.08793 1195.71729
## movement_M10_vmu:phaseweek 2                      0.11231    0.08711 1196.12480
## movement_M10_vmu:phaseweek 3                      0.13962    0.09317 1197.22463
## movement_M10_vmu:phaseweek 4                      0.06966    0.09559 1194.70554
## Allocationlithium:movement_M10_vmu:phaseweek 1   -0.28495    0.10686 1194.77906
## Allocationlithium:movement_M10_vmu:phaseweek 2   -0.22793    0.10675 1195.21961
## Allocationlithium:movement_M10_vmu:phaseweek 3   -0.43293    0.11161 1196.60859
## Allocationlithium:movement_M10_vmu:phaseweek 4   -0.33428    0.11352 1195.18760
##                                                t value Pr(>|t|)    
## (Intercept)                                     -5.048 3.41e-06 ***
## Allocationlithium                                1.177 0.240346    
## movement_M10_vmu                                 5.518 4.20e-08 ***
## phaseweek 1                                      0.093 0.926234    
## phaseweek 2                                      0.082 0.934577    
## phaseweek 3                                     -0.282 0.777814    
## phaseweek 4                                     -1.253 0.210414    
## Age                                             -1.465 0.153452    
## GenderM                                          0.346 0.731626    
## Allocationlithium:movement_M10_vmu               1.890 0.058965 .  
## Allocationlithium:phaseweek 1                   -2.523 0.011773 *  
## Allocationlithium:phaseweek 2                   -1.973 0.048732 *  
## Allocationlithium:phaseweek 3                   -3.228 0.001279 ** 
## Allocationlithium:phaseweek 4                   -2.297 0.021807 *  
## movement_M10_vmu:phaseweek 1                     0.788 0.430709    
## movement_M10_vmu:phaseweek 2                     1.289 0.197536    
## movement_M10_vmu:phaseweek 3                     1.499 0.134256    
## movement_M10_vmu:phaseweek 4                     0.729 0.466259    
## Allocationlithium:movement_M10_vmu:phaseweek 1  -2.667 0.007765 ** 
## Allocationlithium:movement_M10_vmu:phaseweek 2  -2.135 0.032945 *  
## Allocationlithium:movement_M10_vmu:phaseweek 3  -3.879 0.000111 ***
## Allocationlithium:movement_M10_vmu:phaseweek 4  -2.945 0.003296 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 22 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

pC1 <- plot_model(cB1, type = "pred", terms = c("movement_M10_vmu","Allocation","phase"), show.data = F)+
  scale_color_manual(values=orangeNavyPalette)+
  scale_fill_manual(values=orangeNavyPalette)+
  labs(title = "",
       x = "daytime movement volatility",
       y = "movement onset time volatility")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

pC1$facet$params$nrow=1

cB2 <- lmer(pos_vmu ~ Allocation * movement_M10_vmu * phase + Age + Gender + (1|ID),data = df_hlm)
summary(cB2)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pos_vmu ~ Allocation * movement_M10_vmu * phase + Age + Gender +  
##     (1 | ID)
##    Data: df_hlm
## 
## REML criterion at convergence: 3021.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5610 -0.6593 -0.0498  0.6285  4.0737 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.8128   0.9016  
##  Residual             0.5784   0.7605  
## Number of obs: 1243, groups:  ID, 34
## 
## Fixed effects:
##                                                  Estimate Std. Error         df
## (Intercept)                                    -1.462e+00  7.004e-01  7.951e+01
## Allocationlithium                              -1.389e+00  6.316e-01  3.524e+02
## movement_M10_vmu                                5.542e-01  1.013e-01  1.197e+03
## phaseweek 1                                    -2.307e+00  7.245e-01  1.198e+03
## phaseweek 2                                    -4.173e+00  7.334e-01  1.199e+03
## phaseweek 3                                    -4.082e+00  8.030e-01  1.200e+03
## phaseweek 4                                    -1.310e+00  8.371e-01  1.197e+03
## Age                                             5.872e-03  1.348e-02  2.998e+01
## GenderM                                        -6.601e-01  3.193e-01  3.013e+01
## Allocationlithium:movement_M10_vmu             -3.097e-01  1.282e-01  1.200e+03
## Allocationlithium:phaseweek 1                  -5.764e-02  8.594e-01  1.197e+03
## Allocationlithium:phaseweek 2                   7.993e-01  8.843e-01  1.198e+03
## Allocationlithium:phaseweek 3                   3.074e+00  9.433e-01  1.200e+03
## Allocationlithium:phaseweek 4                   1.813e-01  9.699e-01  1.197e+03
## movement_M10_vmu:phaseweek 1                   -4.397e-01  1.514e-01  1.197e+03
## movement_M10_vmu:phaseweek 2                   -7.595e-01  1.499e-01  1.197e+03
## movement_M10_vmu:phaseweek 3                   -7.462e-01  1.604e-01  1.198e+03
## movement_M10_vmu:phaseweek 4                   -2.459e-01  1.645e-01  1.195e+03
## Allocationlithium:movement_M10_vmu:phaseweek 1 -5.375e-03  1.839e-01  1.195e+03
## Allocationlithium:movement_M10_vmu:phaseweek 2  1.364e-01  1.837e-01  1.196e+03
## Allocationlithium:movement_M10_vmu:phaseweek 3  5.283e-01  1.921e-01  1.198e+03
## Allocationlithium:movement_M10_vmu:phaseweek 4 -5.125e-02  1.954e-01  1.196e+03
##                                                t value Pr(>|t|)    
## (Intercept)                                     -2.088  0.04001 *  
## Allocationlithium                               -2.199  0.02851 *  
## movement_M10_vmu                                 5.472 5.42e-08 ***
## phaseweek 1                                     -3.185  0.00149 ** 
## phaseweek 2                                     -5.689 1.60e-08 ***
## phaseweek 3                                     -5.084 4.29e-07 ***
## phaseweek 4                                     -1.564  0.11798    
## Age                                              0.436  0.66627    
## GenderM                                         -2.068  0.04735 *  
## Allocationlithium:movement_M10_vmu              -2.416  0.01586 *  
## Allocationlithium:phaseweek 1                   -0.067  0.94654    
## Allocationlithium:phaseweek 2                    0.904  0.36624    
## Allocationlithium:phaseweek 3                    3.259  0.00115 ** 
## Allocationlithium:phaseweek 4                    0.187  0.85174    
## movement_M10_vmu:phaseweek 1                    -2.905  0.00374 ** 
## movement_M10_vmu:phaseweek 2                    -5.066 4.71e-07 ***
## movement_M10_vmu:phaseweek 3                    -4.653 3.63e-06 ***
## movement_M10_vmu:phaseweek 4                    -1.494  0.13535    
## Allocationlithium:movement_M10_vmu:phaseweek 1  -0.029  0.97669    
## Allocationlithium:movement_M10_vmu:phaseweek 2   0.743  0.45790    
## Allocationlithium:movement_M10_vmu:phaseweek 3   2.750  0.00605 ** 
## Allocationlithium:movement_M10_vmu:phaseweek 4  -0.262  0.79316    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 22 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

pC2 <- plot_model(cB2, type = "pred", terms = c("movement_M10_vmu","Allocation","phase"), show.data = F)+
  scale_color_manual(values=orangeNavyPalette)+
  scale_fill_manual(values=orangeNavyPalette)+
  labs(title = "",
       x = "daytime movement volatility",
       y = "PA volatility")+
  theme(legend.position = "none")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

pC2$facet$params$nrow=1

cB3 <- lmer(pos_vmu ~ Allocation * time_M10_vmu * phase + Age + Gender + (1|ID),data = df_hlm)
summary(cB3)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pos_vmu ~ Allocation * time_M10_vmu * phase + Age + Gender +  
##     (1 | ID)
##    Data: df_hlm
## 
## REML criterion at convergence: 2973.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4155 -0.6696 -0.0367  0.6575  4.1338 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.7574   0.8703  
##  Residual             0.5569   0.7463  
## Number of obs: 1243, groups:  ID, 34
## 
## Fixed effects:
##                                              Estimate Std. Error         df
## (Intercept)                                -2.295e+00  7.187e-01  1.003e+02
## Allocationlithium                           6.977e-01  7.335e-01  5.828e+02
## time_M10_vmu                                4.161e-01  1.244e-01  1.204e+03
## phaseweek 1                                -3.048e-01  6.428e-01  1.196e+03
## phaseweek 2                                -1.291e+00  6.320e-01  1.195e+03
## phaseweek 3                                -1.762e+00  7.224e-01  1.199e+03
## phaseweek 4                                -1.823e+00  7.705e-01  1.200e+03
## Age                                         9.570e-03  1.304e-02  3.023e+01
## GenderM                                    -6.622e-01  3.080e-01  3.004e+01
## Allocationlithium:time_M10_vmu              1.716e-01  1.662e-01  1.202e+03
## Allocationlithium:phaseweek 1              -3.211e+00  9.105e-01  1.199e+03
## Allocationlithium:phaseweek 2              -3.434e+00  8.991e-01  1.199e+03
## Allocationlithium:phaseweek 3              -2.523e+00  9.219e-01  1.199e+03
## Allocationlithium:phaseweek 4              -3.797e+00  9.706e-01  1.200e+03
## time_M10_vmu:phaseweek 1                   -2.202e-03  1.506e-01  1.194e+03
## time_M10_vmu:phaseweek 2                   -1.855e-01  1.449e-01  1.194e+03
## time_M10_vmu:phaseweek 3                   -3.001e-01  1.558e-01  1.196e+03
## time_M10_vmu:phaseweek 4                   -3.314e-01  1.619e-01  1.196e+03
## Allocationlithium:time_M10_vmu:phaseweek 1 -7.456e-01  2.137e-01  1.198e+03
## Allocationlithium:time_M10_vmu:phaseweek 2 -7.923e-01  2.050e-01  1.197e+03
## Allocationlithium:time_M10_vmu:phaseweek 3 -6.578e-01  2.042e-01  1.196e+03
## Allocationlithium:time_M10_vmu:phaseweek 4 -9.397e-01  2.101e-01  1.197e+03
##                                            t value Pr(>|t|)    
## (Intercept)                                 -3.193 0.001881 ** 
## Allocationlithium                            0.951 0.341889    
## time_M10_vmu                                 3.344 0.000852 ***
## phaseweek 1                                 -0.474 0.635452    
## phaseweek 2                                 -2.042 0.041344 *  
## phaseweek 3                                 -2.438 0.014893 *  
## phaseweek 4                                 -2.366 0.018156 *  
## Age                                          0.734 0.468848    
## GenderM                                     -2.150 0.039726 *  
## Allocationlithium:time_M10_vmu               1.032 0.302095    
## Allocationlithium:phaseweek 1               -3.526 0.000437 ***
## Allocationlithium:phaseweek 2               -3.819 0.000141 ***
## Allocationlithium:phaseweek 3               -2.737 0.006287 ** 
## Allocationlithium:phaseweek 4               -3.912 9.67e-05 ***
## time_M10_vmu:phaseweek 1                    -0.015 0.988340    
## time_M10_vmu:phaseweek 2                    -1.281 0.200573    
## time_M10_vmu:phaseweek 3                    -1.926 0.054343 .  
## time_M10_vmu:phaseweek 4                    -2.047 0.040869 *  
## Allocationlithium:time_M10_vmu:phaseweek 1  -3.490 0.000501 ***
## Allocationlithium:time_M10_vmu:phaseweek 2  -3.865 0.000117 ***
## Allocationlithium:time_M10_vmu:phaseweek 3  -3.221 0.001311 ** 
## Allocationlithium:time_M10_vmu:phaseweek 4  -4.474 8.43e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 22 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

pC3 <- plot_model(cB3, type = "pred", terms = c("time_M10_vmu","Allocation","phase"), show.data = F)+
  scale_color_manual(values=orangeNavyPalette)+
  scale_fill_manual(values=orangeNavyPalette)+
  labs(title = "",
       x = "daytime movement onset time volatility",
       y = "PA volatility")+
  theme(legend.position = "none")

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

pC3$facet$params$nrow=1

pC1+pC2+pC3+
  plot_layout(nrow = 3)+
  plot_annotation(tag_levels = c("A"))

ggsave("../figs/FIG5_change_slope.png",width = 7 ,height=8,dpi = 300)